Changeset 638

Timestamp:

09/27/09 00:58:04 (15 years ago)

Author:

smidl

Message:

redesign of PF and MPF to be more flexible and share more code

Location:

library/bdm/estim

Files:

• particles.cpp (modified) (2 diffs)
• particles.h (modified) (7 diffs)

library/bdm/estim/particles.cpp

@@ -5 +5 @@
 using std::endl;
 
-void PF::bayes ( const vec &dt ) {
-        int i;
-        vec lls ( n );
-        ivec ind;
-        double mlls = -std::numeric_limits<double>::infinity(), sum = 0.0;
+void PF::bayes_gensmp(){
+        #pragma parallel for
+        for (int i = 0; i < n; i++ ) {
+                _samples ( i ) = par->samplecond ( _samples ( i ) );
+                lls(i) = 0;
+        }
+}
 
-        for ( i = 0; i < n; i++ ) {
-                //generate new samples from paramater evolution model;
-                vec old_smp=_samples ( i );
-                _samples ( i ) = par->samplecond ( old_smp );
-                lls ( i ) = par->evallogcond ( _samples ( i ), old_smp );
-                lls ( i ) *= obs->evallogcond ( dt, _samples ( i ) );
-
-                if ( lls ( i ) > mlls ) mlls = lls ( i ); //find maximum
-        }
-
+void PF::bayes_weights(){
+// 
+        double mlls=max(lls);
         // compute weights
-        for ( i = 0; i < n; i++ ) {
+        for (int  i = 0; i < n; i++ ) {
                 _w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood
         }
 
         //renormalize
+        double sw=sum(_w);
+        if (!isfinite(sw)) {
+                for (int i=0;i<n;i++){
+                        if (!isfinite(_w(i))) {_w(i)=0;}
+                }
+                sw = sum(_w);
+                if (!isfinite(sw)) {
+                        bdm_error("Particle filter is lost; no particle is good enough.");
+                }
+        }
+        _w /= sw;
+}
+
+void PF::bayes ( const vec &dt ) {
+        int i;
+        // generate samples - time step
+        bayes_gensmp();
+        // weight them - data step
         for ( i = 0; i < n; i++ ) {
-                sum += _w ( i );
-        };
+                lls ( i ) += obs->evallogcond ( dt, _samples ( i ) ); //+= because lls may have something from gensmp!
+        }
 
-        _w ( i ) /= sum; //?
+        bayes_weights();
 
-        ind = est.resample ( resmethod );
+        if (do_resampling()) {
+                est.resample ( resmethod);
+        }
 
 }
@@ -45 +60 @@
 // }
 
+void MPF::bayes ( const vec &dt ) {     
+        // follows PF::bayes in most places!!!
+        
+        int i;
+        int n=pf->__w().length();
+        vec &lls = pf->_lls();
+        
+        // generate samples - time step
+        pf->bayes_gensmp();
+        // weight them - data step
+        #pragma parallel for
+        for ( i = 0; i < n; i++ ) {
+                BMs(i) -> condition(pf->posterior()._sample(i));
+                BMs(i) -> bayes(dt);
+                lls ( i ) += BMs(i)->_ll();
+        }
+
+        pf->bayes_weights();
+        
+        ivec ind;
+        if (pf->do_resampling()) {
+                pf->resample ( ind);
+                
+                #pragma omp parallel for
+                for ( i = 0; i < n; i++ ) {
+                        if ( ind ( i ) != i ) {//replace the current Bm by a new one
+                                delete BMs ( i );
+                                BMs ( i ) = BMs ( ind ( i ) )->_copy_(); //copy constructor
+                        }
+                };
+        }
+};
+
 
 }

library/bdm/estim/particles.h

@@ -36 +36 @@
         Array<vec> &_samples;
         //! Parameter evolution model
-        mpdf *par;
+        shared_ptr<mpdf> par;
         //! Observation model
-        mpdf *obs;
-
+        shared_ptr<mpdf> obs;
+        //! internal structure storing loglikelihood of predictions
+        vec lls;
+        
         //! which resampling method will be used
         RESAMPLING_METHOD resmethod;
-
+        //! resampling threshold; in this case its meaning is minimum ratio of active particles
+        //! For example, for 0.5 resampling is performed when the numebr of active aprticles drops belo 50%.
+        double res_threshold;
+        
         //! \name Options
         //!@{
@@ -58 +63 @@
                 LIDs.set_size ( 5 );
         };
-        /*      PF ( mpdf *par0, mpdf *obs0, epdf *epdf0, int n0 ) :
-                                est ( ),_w ( est._w() ),_samples ( est._samples() ),opt_L_smp(false), opt_L_wei(false)
-                { set_parameters ( par0,obs0,n0 ); set_statistics ( ones ( n0 ),epdf0 ); };*/
-        void set_parameters ( mpdf *par0, mpdf *obs0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
+        
+        void set_parameters (int n0, double res_th0=0.5, RESAMPLING_METHOD rm = SYSTEMATIC ) {
+                n = n0;
+                res_threshold = res_th0;
+                resmethod = rm;
+        };
+        void set_model ( shared_ptr<mpdf> par0, shared_ptr<mpdf> obs0) {
                 par = par0;
                 obs = obs0;
-                n = n0;
-                resmethod = rm;
+                // set values for posterior
+                est.set_rv(par->_rv());
         };
         void set_statistics ( const vec w0, const epdf &epdf0 ) {
                 est.set_statistics ( w0, epdf0 );
         };
+        void set_statistics ( const eEmp &epdf0 ) {
+                bdm_assert_debug(epdf0._rv().equal(par->_rv()),"Incompatibel input");
+                est=epdf0;
+        };
         //!@}
         //! Set posterior density by sampling from epdf0
-//      void set_est ( const epdf &epdf0 );
+        //! Extends original BM::set_options by two more options:
+        //! \li logweights - meaning that all weightes will be logged
+        //! \li logsamples - all samples will be also logged
         void set_options ( const string &opt ) {
                 BM::set_options ( opt );
@@ -78 +92 @@
                 opt_L_smp = ( opt.find ( "logsamples" ) != string::npos );
         }
+        //! bayes I - generate samples and add their weights to lls
+        virtual void bayes_gensmp();
+        //! bayes II - compute weights of the 
+        virtual void bayes_weights();
+        //! important part of particle filtering - decide if it is time to perform resampling
+        virtual bool do_resampling(){   
+                double eff = 1.0 / ( _w * _w );
+                return eff < ( res_threshold*n );
+        }
         void bayes ( const vec &dt );
         //!access function
-        vec* __w() {
-                return &_w;
+        vec& __w() { return _w; }
+        //!access function
+        vec& _lls() { return lls; }
+        RESAMPLING_METHOD _resmethod() const { return resmethod; }
+        //!access function
+        const eEmp& posterior() const {return est;}
+        
+        /*! configuration structure for basic PF
+        \code
+        parameter_pdf   = mpdf_class;         // parameter evolution pdf
+        observation_pdf = mpdf_class;         // observation pdf
+        prior           = epdf_class;         // prior probability density
+        --- optional ---
+        n               = 10;                 // number of particles
+        resmethod       = 'systematic', or 'multinomial', or 'stratified'
+                                                                                  // resampling method
+        res_threshold   = 0.5;                // resample when active particles drop below 50%
+        \endcode
+        */
+        void from_setting(const Setting &set){
+                par = UI::build<mpdf>(set,"parameter_pdf",UI::compulsory);
+                obs = UI::build<mpdf>(set,"observation_pdf",UI::compulsory);
+                
+                prior_from_set(set);
+                resmethod_from_set(set);
+                // set resampling method
+                //set drv
+                //find potential input - what remains in rvc when we subtract rv
+                RV u = par->_rvc().remove_time().subt( par->_rv() ); 
+                //find potential input - what remains in rvc when we subtract x_t
+                RV obs_u = obs->_rvc().remove_time().subt( par->_rv() ); 
+                
+                u.add(obs_u); // join both u, and check if they do not overlap
+                
+                set_drv(concat(obs->_rv(),u) );
+        }
+        //! auxiliary function reading parameter 'resmethod' from configuration file
+        void resmethod_from_set(const Setting &set){
+                string resmeth;
+                if (UI::get(resmeth,set,"resmethod",UI::optional)){
+                        if (resmeth=="systematic") {
+                                resmethod= SYSTEMATIC;
+                        } else  {
+                                if (resmeth=="multinomial"){
+                                        resmethod=MULTINOMIAL;
+                                } else {
+                                        if (resmeth=="stratified"){
+                                                resmethod= STRATIFIED;
+                                        } else {
+                                                bdm_error("Unknown resampling method");
+                                        }
+                                }
+                        }
+                } else {
+                        resmethod=SYSTEMATIC;
+                };
+                if(!UI::get(res_threshold, set, "res_threshold", UI::optional)){
+                        res_threshold=0.5;
+                }
+        }
+        //! load prior information from set and set internal structures accordingly
+        void prior_from_set(const Setting & set){
+                shared_ptr<epdf> pri = UI::build<epdf>(set,"prior",UI::compulsory);
+                
+                eEmp *test_emp=dynamic_cast<eEmp*>(&(*pri));
+                if (test_emp) { // given pdf is sampled
+                        est=*test_emp;
+                } else {
+                        int n;
+                        if (!UI::get(n,set,"n",UI::optional)){n=10;}
+                        // sample from prior
+                        set_statistics(ones(n)/n, *pri);
+                }
+                //validate();
+        }
+        
+        void validate(){
+                n=_w.length();
+                lls=zeros(n);
+                if (par->_rv()._dsize()>0) {
+                        bdm_assert(par->_rv()._dsize()==est.dimension(),"Mismatch of RV and dimension of posterior" );
+                }
+        }
+        //! resample posterior density (from outside - see MPF)
+        void resample(ivec &ind){
+                est.resample(ind,resmethod);
         }
 };
+UIREGISTER(PF);
 
 /*!
 \brief Marginalized Particle filter
 
-Trivial version: proposal = parameter evolution, observation model is not used. (it is assumed to be part of BM).
+A composition of particle filter with exact (or almost exact) bayesian models (BMs).
+The Bayesian models provide marginalized predictive density. Internaly this is achieved by virtual class MPFmpdf.
 */
 
-template<class BM_T>
-class MPF : public PF {
-        Array<BM_T*> BMs;
+class MPF : public BM  {
+    shared_ptr<PF> pf;
+        Array<BM*> BMs;
 
         //! internal class for MPDF providing composition of eEmp with external components
 
         class mpfepdf : public epdf  {
-        protected:
-                eEmp &E;
-                vec &_w;
-                Array<const epdf*> Coms;
+                shared_ptr<PF> &pf;
+                Array<BM*> &BMs;
         public:
-                mpfepdf ( eEmp &E0 ) :
-                                epdf ( ), E ( E0 ),  _w ( E._w() ),
-                                Coms ( _w.length() ) {
-                };
-                //! read statistics from MPF
-                void read_statistics ( Array<BM_T*> &A ) {
-                        dim = E.dimension() + A ( 0 )->posterior().dimension();
-                        for ( int i = 0; i < _w.length() ; i++ ) {
-                                Coms ( i ) = &(A ( i )->posterior());
+                mpfepdf (shared_ptr<PF> &pf0, Array<BM*> &BMs0): epdf(), pf(pf0), BMs(BMs0) { };
+                //! a variant of set parameters - this time, parameters are read from BMs and pf
+                void read_parameters(){
+                        rv = concat(pf->posterior()._rv(), BMs(0)->posterior()._rv());
+                        dim = pf->posterior().dimension() + BMs(0)->posterior().dimension();
+                        bdm_assert_debug(dim == rv._dsize(), "Wrong name ");
+                }
+                vec mean() const {
+                        const vec &w = pf->posterior()._w();
+                        vec pom = zeros ( BMs(0)->posterior ().dimension() );
+                        //compute mean of BMs
+                        for ( int i = 0; i < w.length(); i++ ) {
+                                pom += BMs ( i )->posterior().mean() * w ( i );
                         }
-                }
-                //! needed in resampling
-                void set_elements ( int &i, double wi, const epdf* ep ) {
-                        _w ( i ) = wi;
-                        Coms ( i ) = ep;
-                };
-
-                void set_parameters ( int n ) {
-                        E.set_parameters ( n, false );
-                        Coms.set_length ( n );
-                }
-                vec mean() const {
-                        // ugly
-                        vec pom = zeros ( Coms ( 0 )->dimension() );
-                        for ( int i = 0; i < _w.length(); i++ ) {
-                                pom += Coms ( i )->mean() * _w ( i );
+                        return concat ( pf->posterior().mean(), pom );
+                }
+                vec variance() const {
+                        const vec &w = pf->posterior()._w();
+                        
+                        vec pom = zeros ( BMs(0)->posterior ().dimension() );
+                        vec pom2 = zeros ( BMs(0)->posterior ().dimension() );
+                        vec mea;
+                        
+                        for ( int i = 0; i < w.length(); i++ ) {
+                                // save current mean 
+                                mea = BMs ( i )->posterior().mean();
+                                pom += mea * w ( i );
+                                //compute variance
+                                pom2 += ( BMs ( i )->posterior().variance() + pow ( mea, 2 ) ) * w ( i );
                         }
-                        return concat ( E.mean(), pom );
-                }
-                vec variance() const {
-                        // ugly
-                        vec pom = zeros ( Coms ( 0 )->dimension() );
-                        vec pom2 = zeros ( Coms ( 0 )->dimension() );
-                        for ( int i = 0; i < _w.length(); i++ ) {
-                                pom += Coms ( i )->mean() * _w ( i );
-                                pom2 += ( Coms ( i )->variance() + pow ( Coms ( i )->mean(), 2 ) ) * _w ( i );
-                        }
-                        return concat ( E.variance(), pom2 - pow ( pom, 2 ) );
-                }
+                        return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) );
+                }
+                
                 void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const {
                         //bounds on particles
                         vec lbp;
                         vec ubp;
-                        E.qbounds ( lbp, ubp );
+                        pf->posterior().qbounds ( lbp, ubp );
 
                         //bounds on Components
-                        int dimC = Coms ( 0 )->dimension();
+                        int dimC = BMs ( 0 )->posterior().dimension();
                         int j;
                         // temporary
@@ -160 +261 @@
                         Ubc = -std::numeric_limits<double>::infinity();
 
-                        for ( int i = 0; i < _w.length(); i++ ) {
+                        for ( int i = 0; i < BMs.length(); i++ ) {
                                 // check Coms
-                                Coms ( i )->qbounds ( lbc, ubc );
+                                BMs ( i )->posterior().qbounds ( lbc, ubc );
+                                //save either minima or maxima
                                 for ( j = 0; j < dimC; j++ ) {
                                         if ( lbc ( j ) < Lbc ( j ) ) {
@@ -195 +297 @@
 public:
         //! Default constructor.
-        MPF () : PF (), jest ( est ) {};
-        void set_parameters ( mpdf *par0, mpdf *obs0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
-                PF::set_parameters ( par0, obs0, n0, rm );
-                jest.set_parameters ( n0 );//duplication of rm
+        MPF () :  jest (pf,BMs) {};
+        void set_parameters ( shared_ptr<mpdf> par0, shared_ptr<mpdf> obs0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
+                pf->set_model ( par0, obs0); 
+                pf->set_parameters(n0, rm );
                 BMs.set_length ( n0 );
         }
-        void set_statistics ( const epdf &epdf0, const BM_T* BMcond0 ) {
-
-                PF::set_statistics ( ones ( n ) / n, epdf0 );
+        void set_BM ( const BM &BMcond0 ) {
+
+                int n=pf->__w().length();
+                BMs.set_length(n);
                 // copy
+                //BMcond0 .condition ( pf->posterior()._sample ( 0 ) );
                 for ( int i = 0; i < n; i++ ) {
-                        BMs ( i ) = new BM_T ( *BMcond0 );
-                        BMs ( i )->condition ( _samples ( i ) );
-                }
-
-                jest.read_statistics ( BMs );
-                //options
+                        BMs ( i ) = BMcond0._copy_();
+                        BMs ( i )->condition ( pf->posterior()._sample ( i ) );
+                }
         };
 
@@ -218 +319 @@
                 return jest;
         }
-        //! Set postrior of \c rvc to samples from epdf0. Statistics of BMs are not re-computed! Use only for initialization!
-        /*      void set_est ( const epdf& epdf0 ) {
-                        PF::set_est ( epdf0 );  // sample params in condition
-                        // copy conditions to BMs
-
-                        for ( int i=0;i<n;i++ ) {BMs(i)->condition ( _samples ( i ) );}
-                }*/
+        //! Extends options understood by BM::set_options by option 
+        //! \li logmeans - meaning 
         void set_options ( const string &opt ) {
-                PF::set_options ( opt );
+                BM::set_options(opt);
                 opt_L_mea = ( opt.find ( "logmeans" ) != string::npos );
         }
@@ -234 +330 @@
                 return BMs ( i );
         }
+        
+        /*! configuration structure for basic PF
+        \code
+        BM              = BM_class;           // Bayesian filtr for analytical part of the model
+        parameter_pdf   = mpdf_class;         // transitional pdf for non-parametric part of the model
+        prior           = epdf_class;         // prior probability density
+        --- optional ---
+        n               = 10;                 // number of particles
+        resmethod       = 'systematic', or 'multinomial', or 'stratified'
+                                                                                  // resampling method
+        \endcode
+        */      
+        void from_setting(const Setting &set){
+                shared_ptr<mpdf> par = UI::build<mpdf>(set,"parameter_pdf",UI::compulsory);
+                shared_ptr<mpdf> obs= new mpdf(); // not used!!
+
+                pf = new PF;
+                // rpior must be set before BM
+                pf->prior_from_set(set);
+                pf->resmethod_from_set(set);
+                pf->set_model(par,obs);
+                
+                shared_ptr<BM> BM0 =UI::build<BM>(set,"BM",UI::compulsory);
+                set_BM(*BM0);
+                
+                string opt;
+                if (UI::get(opt,set,"options",UI::optional)){
+                        set_options(opt);
+                }
+                //set drv
+                //find potential input - what remains in rvc when we subtract rv
+                RV u = par->_rvc().remove_time().subt( par->_rv() );            
+                set_drv(concat(BM0->_drv(),u) );
+                validate();
+        }
+        void validate(){
+                try{
+                pf->validate();
+                } catch (std::exception &e){
+                        throw UIException("Error in PF part of MPF:");
+                }
+                jest.read_parameters();
+        }
+        
 };
-
-template<class BM_T>
-void MPF<BM_T>::bayes ( const vec &dt ) {
-        int i;
-        vec lls ( n );
-        vec llsP ( n );
-        ivec ind;
-        double mlls = -std::numeric_limits<double>::infinity();
-
-#pragma omp parallel for
-        for ( i = 0; i < n; i++ ) {
-                //generate new samples from paramater evolution model;
-                vec old_smp=_samples ( i );
-                _samples ( i ) = par->samplecond ( old_smp );
-                llsP ( i ) = par->evallogcond ( _samples ( i ), old_smp );
-                BMs ( i )->condition ( _samples ( i ) );
-                BMs ( i )->bayes ( dt );
-                lls ( i ) = BMs ( i )->_ll(); // lls above is also in proposal her must be lls(i) =, not +=!!
-                if ( lls ( i ) > mlls ) mlls = lls ( i ); //find maximum likelihood (for numerical stability)
-        }
-
-        double sum_w = 0.0;
-        // compute weights
-#pragma omp parallel for
-        for ( i = 0; i < n; i++ ) {
-                _w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood
-                sum_w += _w ( i );
-        }
-
-        if ( sum_w  > 0.0 ) {
-                _w /= sum_w; //?
-        } else {
-                cout << "sum(w)==0" << endl;
-        }
-
-
-        double eff = 1.0 / ( _w * _w );
-        if ( eff < ( 0.3*n ) ) {
-                ind = est.resample ( resmethod );
-                // Resample Bms!
-
-#pragma omp parallel for
-                for ( i = 0; i < n; i++ ) {
-                        if ( ind ( i ) != i ) {//replace the current Bm by a new one
-                                //fixme this would require new assignment operator
-                                // *Bms[i] = *Bms[ind ( i ) ];
-
-                                // poor-man's solution: replicate constructor here
-                                // copied from MPF::MPF
-                                delete BMs ( i );
-                                BMs ( i ) = new BM_T ( *BMs ( ind ( i ) ) ); //copy constructor
-                                const epdf& pom = BMs ( i )->posterior();
-                                jest.set_elements ( i, 1.0 / n, &pom );
-                        }
-                };
-                cout << '.';
-        }
-}
+UIREGISTER(MPF);
 
 }