root/library/bdm/estim/particles.cpp @ 661

#include "particles.h"

namespace bdm {

using std::endl;

void PF::bayes_gensmp(){
        #pragma parallel for
        for (int i = 0; i < n; i++ ) {
                _samples ( i ) = par->samplecond ( _samples ( i ) );
                lls(i) = 0;
        }
}

void PF::bayes_weights(){
// 
        double mlls=max(lls);
        // compute weights
        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++ ) {
                lls ( i ) += obs->evallogcond ( dt, _samples ( i ) ); //+= because lls may have something from gensmp!
        }

        bayes_weights();

        if (do_resampling()) {
                est.resample ( resmethod);
        }

}

// void PF::set_est ( const epdf &epdf0 ) {
//      int i;
//
//      for ( i=0;i<n;i++ ) {
//              _samples ( i ) = epdf0.sample();
//      }
// }

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
                        }
                };
        }
};


}
//MPF::MPF:{}