Changeset 679 for library/bdm/estim/kalman.h

Timestamp:

10/23/09 00:05:25 (15 years ago)

Author:

smidl

Message:

Major changes in BM -- OK is only test suite and tests/tutorial -- the rest is broken!!!

Files:

• library/bdm/estim/kalman.h (modified) (17 diffs)

library/bdm/estim/kalman.h

@@ -34 +34 @@
 {
         protected:
-                //! cache of rv.count()
-                int dimx;
-                //! cache of rvy.count()
-                int dimy;
-                //! cache of rvu.count()
-                int dimu;
                 //! Matrix A
                 mat A;
@@ -53 +47 @@
                 sq_T R;
         public:
-                StateSpace() : dimx (0), dimy (0), dimu (0), A(), B(), C(), D(), Q(), R() {}
+                StateSpace() :  A(), B(), C(), D(), Q(), R() {}
                 //!copy constructor
-                StateSpace(const StateSpace<sq_T> &S0) : dimx (S0.dimx), dimy (S0.dimy), dimu (S0.dimu), A(S0.A), B(S0.B), C(S0.C), D(S0.D), Q(S0.Q), R(S0.R) {}
+                StateSpace(const StateSpace<sq_T> &S0) :  A(S0.A), B(S0.B), C(S0.C), D(S0.D), Q(S0.Q), R(S0.R) {}
                 //! set all matrix parameters
                 void set_parameters (const mat &A0, const  mat &B0, const  mat &C0, const  mat &D0, const  sq_T &Q0, const sq_T &R0);
@@ -83 +77 @@
                 }               
                 //! access function
-                int _dimx(){return dimx;}
-                //! access function
-                int _dimy(){return dimy;}
-                //! access function
-                int _dimu(){return dimu;}
-                //! access function
                 const mat& _A() const {return A;}
                 //! access function
@@ -109 +97 @@
                 //! id of output
                 RV yrv;
-                //! id of input
-                RV urv;
                 //! Kalman gain
                 mat  _K;
                 //!posterior
-                shared_ptr<enorm<sq_T> > est;
+                enorm<sq_T> est;
                 //!marginal on data f(y|y)
                 enorm<sq_T>  fy;
         public:
-                Kalman<sq_T>() : BM(), StateSpace<sq_T>(), yrv(),urv(), _K(),  est(new enorm<sq_T>){}
+                Kalman<sq_T>() : BM(), StateSpace<sq_T>(), yrv(), _K(),  est(){}
                 //! Copy constructor
-                Kalman<sq_T>(const Kalman<sq_T> &K0) : BM(K0), StateSpace<sq_T>(K0), yrv(K0.yrv),urv(K0.urv), _K(K0._K),  est(new enorm<sq_T>(*K0.est)), fy(K0.fy){}
+                Kalman<sq_T>(const Kalman<sq_T> &K0) : BM(K0), StateSpace<sq_T>(K0), yrv(K0.yrv), _K(K0._K),  est(K0.est), fy(K0.fy){}
                 //!set statistics of the posterior
-                void set_statistics (const vec &mu0, const mat &P0) {est->set_parameters (mu0, P0); };
+                void set_statistics (const vec &mu0, const mat &P0) {est.set_parameters (mu0, P0); };
                 //!set statistics of the posterior
-                void set_statistics (const vec &mu0, const sq_T &P0) {est->set_parameters (mu0, P0); };
+                void set_statistics (const vec &mu0, const sq_T &P0) {est.set_parameters (mu0, P0); };
                 //! return correctly typed posterior (covariant return)
-                const enorm<sq_T>& posterior() const {return *est.get();}
-                //! shared posterior
-                shared_ptr<epdf> shared_posterior() {return est;}
+                const enorm<sq_T>& posterior() const {return est;}
                 //! load basic elements of Kalman from structure
                 void from_setting (const Setting &set) {
@@ -139 +123 @@
                         // Initial values
                         UI::get (yrv, set, "yrv", UI::optional);
-                        UI::get (urv, set, "urv", UI::optional);
-                        set_drv(concat(yrv,urv));
+                        UI::get (rvc, set, "urv", UI::optional);
+                        set_yrv(concat(yrv,rvc));
                         
                         validate();
@@ -147 +131 @@
                 void validate() {
                         StateSpace<sq_T>::validate();
-                        bdm_assert(est->dimension(), "Statistics and model parameters mismatch");
+                        dimy = this->C.rows();
+                        dimc = this->B.cols();
+                        set_dim(this->A.rows());
+                        
+                        bdm_assert(est.dimension(), "Statistics and model parameters mismatch");
                 }
 };
@@ -160 +148 @@
                 KalmanFull() :Kalman<fsqmat>(){};
                 //! Here dt = [yt;ut] of appropriate dimensions
-                void bayes (const vec &dt);
+                void bayes (const vec &yt, const vec &cond=empty_vec);
                 BM* _copy_() const {
                         KalmanFull* K = new KalmanFull;
                         K->set_parameters (A, B, C, D, Q, R);
-                        K->set_statistics (est->_mu(), est->_R());
+                        K->set_statistics (est._mu(), est._R());
                         return K;
                 }
@@ -192 +180 @@
                         KalmanCh* K = new KalmanCh;
                         K->set_parameters (A, B, C, D, Q, R);
-                        K->set_statistics (est->_mu(), est->_R());
+                        K->set_statistics (est._mu(), est._R());
                         return K;
                 }
@@ -213 +201 @@
                 Thus this object evaluates only predictors! Not filtering densities.
                 */
-                void bayes (const vec &dt);
+                void bayes (const vec &yt, const vec &cond=empty_vec);
 
                 void from_setting(const Setting &set){
                         Kalman<chmat>::from_setting(set);
+                        validate();
+                }
+                void validate() {
+                        Kalman<chmat>::validate();
                         initialize();
                 }
@@ -244 +236 @@
 
                 //! Here dt = [yt;ut] of appropriate dimensions
-                void bayes (const vec &dt);
+                void bayes (const vec &yt, const vec &cond=empty_vec);
                 //! set estimates
                 void set_statistics (const vec &mu0, const mat &P0) {
-                        est->set_parameters (mu0, P0);
+                        est.set_parameters (mu0, P0);
                 };
                 //! access function
                 const mat _R() {
-                        return est->_R().to_mat();
+                        return est._R().to_mat();
                 }
                 void from_setting (const Setting &set) {
@@ -280 +272 @@
                         //connect
                         shared_ptr<RV> drv = UI::build<RV> ( set, "drv", UI::compulsory );
-                        set_drv ( *drv );
+                        set_yrv ( *drv );
                         shared_ptr<RV> rv = UI::build<RV> ( set, "rv", UI::compulsory );
                         set_rv ( *rv );
@@ -337 +329 @@
                         EKFCh* E = new EKFCh;
                         E->set_parameters (pfxu, phxu, Q, R);
-                        E->set_statistics (est->_mu(), est->_R());
+                        E->set_statistics (est._mu(), est._R());
                         return E;
                 }
@@ -344 +336 @@
 
                 //! Here dt = [yt;ut] of appropriate dimensions
-                void bayes (const vec &dt);
+                void bayes (const vec &yt, const vec &cond=empty_vec);
 
                 void from_setting (const Setting &set);
@@ -389 +381 @@
                         est.set_parameters (A (0)->posterior().mean(), A (0)->posterior()._R());
                 }
-                void bayes (const vec &dt) {
+                void bayes (const vec &yt, const vec &cond=empty_vec) {
                         int n = Models.length();
                         int i;
                         for (i = 0; i < n; i++) {
-                                Models (i)->bayes (dt);
+                                Models (i)->bayes (yt);
                                 _lls (i) = Models (i)->_ll();
                         }
@@ -475 +467 @@
                                 Crv.add(urv.copy_t(t));
                         }
-                        
-                        this->dimx = xrv._dsize();
-                        this->dimy = yrv._dsize();
-                        this->dimu = urv._dsize();
-                        
+                                                
                         // get mapp
                         th2A.set_connection(xrv, ml._rvc());
@@ -486 +474 @@
 
                         //set matrix sizes
-                        this->A=zeros(dimx,dimx);
-                        for (int j=1; j<dimx; j++){A(j,j-1)=1.0;} // off diagonal
-                                this->B=zeros(dimx,1);
+                        this->A=zeros(xrv._dsize(),xrv._dsize());
+                        for (int j=1; j<xrv._dsize(); j++){A(j,j-1)=1.0;} // off diagonal
+                                this->B=zeros(xrv._dsize(),1);
                                 this->B(0) = 1.0;
-                                this->C=zeros(1,dimx);
+                                this->C=zeros(1,xrv._dsize());
                                 this->D=zeros(1,urv._dsize());
-                                this->Q = zeros(dimx,dimx);
+                                this->Q = zeros(xrv._dsize(),xrv._dsize());
                         // R is set by update
                         
@@ -538 +526 @@
 template<class sq_T>
 void StateSpace<sq_T>::validate(){
-        dimx = A.rows();
-        dimu = B.cols();
-        dimy = C.rows();
-        bdm_assert (A.cols() == dimx, "KalmanFull: A is not square");
-        bdm_assert (B.rows() == dimx, "KalmanFull: B is not compatible");
-        bdm_assert (C.cols() == dimx, "KalmanFull: C is not square");
-        bdm_assert ( (D.rows() == dimy) && (D.cols() == dimu), "KalmanFull: D is not compatible");
-        bdm_assert ( (Q.cols() == dimx) && (Q.rows() == dimx), "KalmanFull: Q is not compatible");
-        bdm_assert ( (R.cols() == dimy) && (R.rows() == dimy), "KalmanFull: R is not compatible");
+        bdm_assert (A.cols() == A.rows(), "KalmanFull: A is not square");
+        bdm_assert (B.rows() == A.rows(), "KalmanFull: B is not compatible");
+        bdm_assert (C.cols() == A.rows(), "KalmanFull: C is not compatible");
+        bdm_assert ( (D.rows() == C.rows()) && (D.cols() == B.cols()), "KalmanFull: D is not compatible");
+        bdm_assert ( (Q.cols() == A.rows()) && (Q.rows() == A.rows()), "KalmanFull: Q is not compatible");
+        bdm_assert ( (R.cols() == C.rows()) && (R.rows() == C.rows()), "KalmanFull: R is not compatible");
 }