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

#include "particles.h"

namespace bdm {

using std::endl;

void PF::bayes_gensmp ( const vec &ut ) {
        if ( ut.length() > 0 ) {
                vec cond ( par->dimensionc() );
                cond.set_subvector ( par->dimension(), ut );
#pragma parallel for
                for ( int i = 0; i < n; i++ ) {
                        cond.set_subvector ( 0, _samples ( i ) );
                        _samples ( i ) = par->samplecond ( cond );
                        lls ( i ) = 0;
                }
        } else {
#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 ( !std::isfinite ( sw ) ) {
                for ( int i = 0; i < n; i++ ) {
                        if ( !std::isfinite ( _w ( i ) ) ) {
                                _w ( i ) = 0;
                        }
                }
                sw = sum ( _w );
                if ( !std::isfinite ( sw ) || sw == 0.0 ) {
                        bdm_error ( "Particle filter is lost; no particle is good enough." );
                }
        }
        _w /= sw;
}

void PF::bayes ( const vec &yt, const vec &cond ) {
        const vec &ut = cond; //todo

        int i;
        // generate samples - time step
        bayes_gensmp ( ut );
        // weight them - data step
        for ( i = 0; i < n; i++ ) {
                lls ( i ) += obs->evallogcond ( yt, _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 &yt, const vec &cond ) {
        // follows PF::bayes in most places!!!
        int i;
        int n = pf->__w().length();
        vec &lls = pf->_lls();
        Array<vec> &samples = pf->__samples();

        // generate samples - time step
        pf->bayes_gensmp ( vec ( 0 ) );
        // weight them - data step
#pragma parallel for
        for ( i = 0; i < n; i++ ) {
                vec bm_cond ( BMs ( i )->dimensionc() );
                this2bm.fill_cond ( yt, cond, bm_cond );
                pf2bm.filldown ( samples ( i ), bm_cond );
                BMs ( i ) -> bayes ( this2bm.pushdown ( yt ), bm_cond );
                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:{}