Changeset 477 for library/bdm/estim/particles.h

Timestamp:

08/05/09 14:40:03 (15 years ago)

Author:

mido

Message:

panove, vite, jak jsem peclivej na upravu kodu.. snad se vam bude libit:) konfigurace je v souboru /system/astylerc

Files:

• library/bdm/estim/particles.h (modified) (12 diffs)

library/bdm/estim/particles.h

@@ -55 +55 @@
         //! \name Constructors
         //!@{
-        PF ( ) :est(), _w ( est._w() ),_samples ( est._samples() ), opt_L_smp ( false ), opt_L_wei ( false ) {LIDs.set_size ( 5 );};
+        PF ( ) : est(), _w ( est._w() ), _samples ( est._samples() ), opt_L_smp ( false ), opt_L_wei ( false ) {
+                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 )
-        { par = par0; obs=obs0; n=n0; resmethod= rm;};
-        void set_statistics ( const vec w0, epdf *epdf0 ) {est.set_statistics ( w0,epdf0 );};
+        void set_parameters ( mpdf *par0, mpdf *obs0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
+                par = par0;
+                obs = obs0;
+                n = n0;
+                resmethod = rm;
+        };
+        void set_statistics ( const vec w0, epdf *epdf0 ) {
+                est.set_statistics ( w0, epdf0 );
+        };
         //!@}
         //! Set posterior density by sampling from epdf0
 //      void set_est ( const epdf &epdf0 );
         void set_options ( const string &opt ) {
-                BM::set_options(opt);
-                opt_L_wei= ( opt.find ( "logweights" ) !=string::npos );
-                opt_L_smp= ( opt.find ( "logsamples" ) !=string::npos );
+                BM::set_options ( opt );
+                opt_L_wei = ( opt.find ( "logweights" ) != string::npos );
+                opt_L_smp = ( opt.find ( "logsamples" ) != string::npos );
         }
         void bayes ( const vec &dt );
         //!access function
-        vec* __w() {return &_w;}
+        vec* __w() {
+                return &_w;
+        }
 };
 
@@ -87 +97 @@
         //! internal class for MPDF providing composition of eEmp with external components
 
-class mpfepdf : public epdf  {
+        class mpfepdf : public epdf  {
         protected:
                 eEmp &E;
@@ -99 +109 @@
                 //! 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 )->_e();}
+                        dim = E.dimension() + A ( 0 )->posterior().dimension();
+                        for ( int i = 0; i < _w.length() ; i++ ) {
+                                Coms ( i ) = A ( i )->_e();
+                        }
                 }
                 //! needed in resampling
-                void set_elements ( int &i, double wi, const epdf* ep )
-                {_w ( i ) =wi; Coms ( i ) =ep;};
+                void set_elements ( int &i, double wi, const epdf* ep ) {
+                        _w ( i ) = wi;
+                        Coms ( i ) = ep;
+                };
 
                 void set_parameters ( int n ) {
@@ -112 +126 @@
                 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 ( E.mean(),pom );
+                        vec pom = zeros ( Coms ( 0 )->dimension() );
+                        for ( int i = 0; i < _w.length(); i++ ) {
+                                pom += Coms ( i )->mean() * _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++ ) {
+                        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 ) );
-                }
-                void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const {
+                                pom2 += ( Coms ( i )->variance() + pow ( Coms ( i )->mean(), 2 ) ) * _w ( i );
+                        }
+                        return concat ( E.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 );
+                        E.qbounds ( lbp, ubp );
 
                         //bounds on Components
-                        int dimC=Coms ( 0 )->dimension();
+                        int dimC = Coms ( 0 )->dimension();
                         int j;
                         // temporary
-                        vec lbc(dimC);
-                        vec ubc(dimC);
+                        vec lbc ( dimC );
+                        vec ubc ( dimC );
                         // minima and maxima
-                        vec Lbc(dimC);
-                        vec Ubc(dimC);
+                        vec Lbc ( dimC );
+                        vec Ubc ( dimC );
                         Lbc = std::numeric_limits<double>::infinity();
                         Ubc = -std::numeric_limits<double>::infinity();
 
-                        for ( int i=0;i<_w.length();i++ ) {
+                        for ( int i = 0; i < _w.length(); i++ ) {
                                 // check Coms
-                                Coms ( i )->qbounds ( lbc,ubc );
-                                for ( j=0;j<dimC; j++ ) {
-                                        if ( lbc ( j ) <Lbc ( j ) ) {Lbc ( j ) =lbc ( j );}
-                                        if ( ubc ( j ) >Ubc ( j ) ) {Ubc ( j ) =ubc ( j );}
+                                Coms ( i )->qbounds ( lbc, ubc );
+                                for ( j = 0; j < dimC; j++ ) {
+                                        if ( lbc ( j ) < Lbc ( j ) ) {
+                                                Lbc ( j ) = lbc ( j );
+                                        }
+                                        if ( ubc ( j ) > Ubc ( j ) ) {
+                                                Ubc ( j ) = ubc ( j );
+                                        }
                                 }
                         }
-                        lb=concat(lbp,Lbc);
-                        ub=concat(ubp,Ubc);
-                }
-
-                vec sample() const {it_error ( "Not implemented" );return 0;}
-
-                double evallog ( const vec &val ) const {it_error ( "not implemented" ); return 0.0;}
+                        lb = concat ( lbp, Lbc );
+                        ub = concat ( ubp, Ubc );
+                }
+
+                vec sample() const {
+                        it_error ( "Not implemented" );
+                        return 0;
+                }
+
+                double evallog ( const vec &val ) const {
+                        it_error ( "not implemented" );
+                        return 0.0;
+                }
         };
 
@@ -170 +196 @@
         //! Default constructor.
         MPF () : PF (), jest ( est ) {};
-        void set_parameters ( mpdf *par0, mpdf *obs0, int n0, RESAMPLING_METHOD rm=SYSTEMATIC ) {
+        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
@@ -177 +203 @@
         void set_statistics ( epdf *epdf0, const BM_T* BMcond0 ) {
 
-                PF::set_statistics ( ones ( n ) /n, epdf0 );
+                PF::set_statistics ( ones ( n ) / n, epdf0 );
                 // copy
-                for ( int i=0;i<n;i++ ) { BMs ( i ) = new BM_T ( *BMcond0 ); BMs ( i )->condition ( _samples ( i ) );}
+                for ( int i = 0; i < n; i++ ) {
+                        BMs ( i ) = new BM_T ( *BMcond0 );
+                        BMs ( i )->condition ( _samples ( i ) );
+                }
 
                 jest.read_statistics ( BMs );
@@ -186 +215 @@
 
         void bayes ( const vec &dt );
-        const epdf& posterior() const {return jest;}
-        const epdf* _e() const {return &jest;} //Fixme: is it useful?
+        const epdf& posterior() const {
+                return jest;
+        }
+        const epdf* _e() const {
+                return &jest;    //Fixme: is it useful?
+        }
         //! 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 ) {
@@ -197 +230 @@
         void set_options ( const string &opt ) {
                 PF::set_options ( opt );
-                opt_L_mea = ( opt.find ( "logmeans" ) !=string::npos );
+                opt_L_mea = ( opt.find ( "logmeans" ) != string::npos );
         }
 
         //!Access function
-        BM* _BM ( int i ) {return BMs ( i );}
+        BM* _BM ( int i ) {
+                return BMs ( i );
+        }
 };
 
@@ -210 +245 @@
         vec llsP ( n );
         ivec ind;
-        double mlls=-std::numeric_limits<double>::infinity();
+        double mlls = -std::numeric_limits<double>::infinity();
 
 #pragma omp parallel for
-        for ( i=0;i<n;i++ ) {
+        for ( i = 0; i < n; i++ ) {
                 //generate new samples from paramater evolution model;
                 _samples ( i ) = par->samplecond ( _samples ( i ) );
@@ -220 +255 @@
                 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;
+                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++ ) {
+        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 );
+                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 );
@@ -245 +279 @@
 
 #pragma omp parallel for
-                for ( i=0;i<n;i++ ) {
-                        if ( ind ( i ) !=i ) {//replace the current Bm by a new one
+                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 ) ];
@@ -254 +288 @@
                                 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 );
+                                const epdf& pom = BMs ( i )->posterior();
+                                jest.set_elements ( i, 1.0 / n, &pom );
                         }
                 };