Changeset 265

Timestamp:

02/09/09 23:14:58 (15 years ago)

Author:

smidl

Message:

UI in matlab

Files:

• bdm/stat/libBM.cpp (modified) (1 diff)
• bdm/stat/libBM.h (modified) (4 diffs)
• bdm/stat/libDS.cpp (modified) (3 diffs)
• bdm/stat/libDS.h (modified) (1 diff)
• matlab/mex/arx_example.cfg (copied) (copied from tests/UI/UIArxDS_test.cfg) (1 diff)
• matlab/mex/build.sh (added)
• matlab/mex/estimator.cpp (added)
• matlab/mex/estimator_disp.m (added)
• matlab/mex/mexlog.h (added)
• tests/UI/UIArxDS_test.cfg (modified) (2 diffs)
• tests/UI/UIArxDS_test.cpp (modified) (2 diffs)

bdm/stat/libBM.cpp

@@ -52 +52 @@
 RV::RV() : tsize ( 0 ), len ( 0 ), ids ( 0 ),  sizes ( 0 ), times ( 0 ),names ( 0 ) {};
 
-RV::RV(string name, int id){
-        it_assert(id>=RVcounter,"RV::This id is already taken");
-        RVcounter=id;
+RV::RV(string name, int id, int sz, int tm){
+        if (id>RVcounter) {RVcounter=id;};
         Array<string> A(1); A(0)=name;
-        init(vec_1(id),A,vec_1(1),vec_1(0));
+        init(vec_1(id),A,vec_1(sz),vec_1(tm));
 }
 

bdm/stat/libBM.h

@@ -78 +78 @@
                 RV ();
                 //! Constructor of a single RV with given id
-                RV (string name, int id);
+                RV (string name, int id, int sz=1, int tm=0);
 
                 //! Printing output e.g. for debugging.
@@ -490 +490 @@
                         L.logit ( L_ut,tmp.mid ( Drv.count(), Urv.count() ) );
                 }
-        };
-
-        /*! \brief Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.
+                //!access function 
+                virtual RV _drv() const {return concat(Drv,Urv);}
+                //!access function 
+                const RV& _urv() const {return Urv;}
+        };
+
+        /*! \brief Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities.
 
         */
@@ -500 +504 @@
                 //!Random variable of the posterior
                 RV rv;
+                //! Random variable of the data (optional)
+                RV drv;
                 //!Logarithm of marginalized data likelihood.
                 double ll;
@@ -537 +543 @@
                 //!access function
                 const RV& _rv() const {return rv;}
+                //!access function 
+                const RV& _drv() const {return drv;}
+                //!set drv
+                void set_drv(const RV &rv){drv=rv;}
                 //!access function
                 double _ll() const {return ll;}

bdm/stat/libDS.cpp

@@ -44 +44 @@
 void ArxDS::step() {
         double tmp;
+        //shift history
+        H.shift_right ( 0, Drv.count() +Urv.count() );
+
+        H.set_subvector ( Drv.count(), U ); // write U after Drv
+
         //get regressor
         rgr = rgrlnk.get_val ( H );
         // Eval Y
-        Y = model.samplecond ( rgr,tmp );
+        H.set_subvector ( 0, model.samplecond ( rgr,tmp ) );
 
-        //shift history
-        H.shift_right ( 0, Drv.count() +Urv.count() );
-        //fill new data
-        H.set_subvector ( Drv.count(),Y );
-
-        //Leaving U.length() empty slots - these will be filled by write()
 }
 
@@ -60 +59 @@
 RV fullrgr ( const RV &drv0, const RV &urv0, const RV &rrv0 ) {
         RV T ( urv0 );
-        RV pom = concat(drv0, urv0);
+        RV pom = concat ( drv0, urv0 );
         int mint = rrv0.mint();
         for ( int i=0; i>mint; i-- ) {
@@ -69 +68 @@
 }
 
-ArxDS::ArxDS ( RV &drv, RV &urv, RV &rrv ) : 
-                DS(drv,urv), Rrv(rrv), Hrv( fullrgr ( drv,urv,rrv )), //RVs
-                Y(drv.count()), H(Hrv.count()) ,rgr ( rrv.count() ),  //tmp variables
-                rgrlnk (Hrv ,rrv ) ,model ( drv,rrv ) {
+ArxDS::ArxDS ( RV &drv, RV &urv, RV &rrv ) :
+                DS ( drv,urv ), Rrv ( rrv ), Hrv ( concat(drv,fullrgr ( drv,urv,rrv )) ), //RVs
+                H ( Hrv.count() ) ,U ( urv.count() ),rgr ( rrv.count() ),  //tmp variables
+                rgrlnk ( rrv, Hrv ) ,model ( drv,rrv ) {
 }

bdm/stat/libDS.h

@@ -53 +53 @@
         */
         class ArxDS : public DS {
+        protected:
                 //! Rv of the regressor
                 RV Rrv;
                 //! Rv of the history (full regressor)
                 RV Hrv;
-                //! Internal storage of results
-                vec Y;
-                //! History, ordered as \f$[u_t, y_{t-1 }, u_{t-1}, \ldots]\f$ 
+                //! History, ordered as \f$[y_t, u_t, y_{t-1 }, u_{t-1}, \ldots]\f$
                 vec H;
+                //! (future) input
+                vec U;
                 //! temporary variable for regressor
                 vec rgr;
-                //! data link: val = y, cond = u; local = rgr;
+                //! data link: H -> rgr
                 datalink_e2e rgrlnk;
-                //! model of Y - linear Gaussian 
+                //! model of Y - linear Gaussian
                 mlnorm<chmat> model;
+                //! options
+                bool opt_L_theta;
+                //! loggers
+                int L_theta;
+                int L_R;
                 public:
-                void getdata ( vec &dt ){it_assert_debug(dt.length()==Y.length(),"ArxDS"); 
-                        dt=concat(Y,H.left(Urv.count()));};
-                void getdata ( vec &dt, const ivec &indexes ){it_assert_debug(dt.length()==Y.length(),"ArxDS"); dt=Y(indexes);};
-                void write ( vec &ut ){it_assert_debug(ut.length()==Urv.count(),"ArxDS"); H.set_subvector(0,ut);};
-                void write ( vec &ut, const ivec &indexes ){it_assert_debug(ut.length()==Urv.count(),"ArxDS"); set_subvector(H,indexes,ut);};
+                void getdata ( vec &dt ) {
+                        it_assert_debug ( dt.length() ==Drv.count(),"ArxDS" );
+                        dt=H.left ( Urv.count() +Drv.count() );
+                };
+                void getdata ( vec &dt, const ivec &indexes ) {
+                        it_assert_debug ( dt.length() ==indeces.length(),"ArxDS" );
+                        dt=H ( indexes );
+                };
+                void write ( vec &ut ) {
+                        it_assert_debug ( ut.length() ==Urv.count(),"ArxDS" );
+                        U=ut;
+                };
+                void write ( vec &ut, const ivec &indexes ) {
+                        it_assert_debug ( ut.length() ==indeces.length(),"ArxDS" );
+                        set_subvector ( U, indexes,ut );
+                };
                 void step();
                 //!Default constructor
-                ArxDS ( RV &drv, RV &urv, RV &rrv);
+                ArxDS ( RV &drv, RV &urv, RV &rrv );
                 //! Set parameters of the internal model
-                void set_parameters(const mat &Th0, const vec mu0, const chmat &sqR0)
-                {model.set_parameters(Th0, mu0, sqR0); };
+                void set_parameters ( const mat &Th0, const vec mu0, const chmat &sqR0 )
+                { model.set_parameters ( Th0, mu0, sqR0 ); };
+                //! set options from a string
+                void set_options ( const string &s ) {
+                        opt_L_theta= ( s.find ( "L_theta" ) !=string::npos );
+                };
+                virtual void log_add ( logger &L ) {
+                        DS::log_add ( L );
+                        mat &A =model._A();
+                        mat R =model._R();
+                        if ( opt_L_theta ) {L_theta=L.add ( RV("{theta }", vec_1(A.rows() *A.cols())),"t" );}
+                        if ( opt_L_theta ) {L_R=L.add ( RV("{R }", vec_1(R.rows() *R.cols())),"r" );}
+                }
+                virtual void logit ( logger &L ) {
+                        DS::logit ( L );
+                        mat &A =model._A();
+                        mat R =model._R();
+                        if ( opt_L_theta ) {L.logit ( L_theta,vec ( A._data(), A.rows() *A.cols() ) );};
+                        if ( opt_L_theta ) {L.logit ( L_R, vec ( R._data(), R.rows() *R.rows() ) );};
+                }
+
         };
 
+        class ARXDS : public ArxDS {
+        public:
+                ARXDS ( RV &drv, RV &urv, RV &rrv ) : ArxDS ( drv,urv,rrv ) {}
+
+                void getdata ( vec &dt ) {dt=H;}
+                void getdata ( vec &dt, const ivec &indeces ) {dt=H ( indeces );}
+                virtual RV _drv() const {return Hrv;}
+
+        };
 }; //namespace
 

matlab/mex/arx_example.cfg

@@ -18 +18 @@
         r = [0.1, 0.0,
              0.0, 1.0];
+        opt="L_theta";
 };
 
 //store results
 logger = {
-        type= "dirfilelog";
-        dirname = "exp/arx_ui";
-        maxlen = 1000; //
+        type= "mexlog";
+        maxlen = 90; //
+    dirname = "exp/ax";
 };
 
 //estimation
 estimator = {
-    type = "ARX";
-        
+   type = "ARXest";
+        ychns = [1];
+        rgrid =  [ 1,  1,  1,  2];
+        delays = [-1, -2, -3, -1];
+
+        //optional fields
+        dV0 = [1.0, 1e-5, 1e-5, 1e-5, 1e-5]; //default: 1e-3 for y, 1e-5 for rgr
+        nu0 = 8.;      //default: rgrlen + 2
+        frg = 1.;    // forgetting, default frg=1.0
 };
 
 //experiment description
 experiment:{
-        ndat = 9000;
+        ndat = 90;
 };

tests/UI/UIArxDS_test.cfg

@@ -18 +18 @@
         r = [0.1, 0.0,
              0.0, 1.0];
+        opt="L_theta";
 };
 
@@ -29 +30 @@
 //estimation
 estimator = {
-    type = "ARX";
-        
+   type = "ARXest";
+        ychns = [1];
+        rgrid =  [ 1,  1,  1,  2];
+        delays = [-1, -2, -3, -1];
+
+        //optional fields
+        dV0 = [1e-3, 1e-5, 1e-5, 1e-5, 1e-5]; //default: 1e-3 for y, 1e-5 for rgr
+        nu0 = 8.;      //default: rgrlen + 2
+        frg = .990;    // forgetting, default frg=1.0
 };
 

tests/UI/UIArxDS_test.cpp

@@ -1 +1 @@
 /*!
-\file 
+\file
 \brief Test of UI builders for ARX
 
@@ -9 +9 @@
 
 #include <stat/libDS_ui.h>
+#include <estim/arx_ui.h>
 
 using namespace bdm;
 int main() {
-        UIFile F("UIArxDS_test.cfg");
+        UIFile F ( "UIArxDS_test.cfg" );
 
-        logger* L; 
-        UIbuild(F.lookup("logger"),L);
+        logger* L;
+        UIbuild ( F.lookup ( "logger" ),L );
         ArxDS * DS;
-        UIbuild(F.lookup("system"),DS);
+        UIbuild ( F.lookup ( "system" ),DS );
         int Ndat;
-        F.lookupValue("experiment.ndat",Ndat);
+        F.lookupValue ( "experiment.ndat",Ndat );
         // SET SIMULATOR
+        BM* E;
+        UIbuild ( F.lookup ( "estimator" ),E );
 
-        DS->log_add(*L);
+        DS->log_add ( *L );
+        int L_est= L->add ( E->_epdf()._rv(), "est" ); // estimate
+        int L_lb = L->add ( E->_epdf()._rv(), "lb" ); // lower bound
+        int L_ub = L->add ( E->_epdf()._rv(), "ub" ); // upper bound
         L->init();
-        
+
+        vec dt=zeros ( DS->_drv().count() );   //data variable
+        datalink_e2e dl ( E->_drv(),DS->_drv() ); //datalink between a datasource and estimator
+
         for ( int tK=1;tK<Ndat;tK++ ) {
-                DS->step();
-                DS->logit(*L);
+                DS->step();                                                     // simulator step
+                DS->getdata ( dt );                                     // read data
+                E->bayes ( dl.get_val ( dt ) );         // update estimates
+
+                DS->logit ( *L );
+                L->logit ( L_est, E->_epdf().mean() );
+                L->logit ( L_lb,  E->_epdf().mean()-2*sqrt ( E->_epdf().variance() ) );
+                L->logit ( L_ub,  E->_epdf().mean() +2*sqrt ( E->_epdf().variance() ) );
+
                 L->step();
         }
-        
-        L->finalize(); 
-        
+
+        L->finalize();
+
         delete L;
         delete DS;
+        delete E;
         return 0;
 }