Changeset 738

Timestamp:

11/25/09 12:46:08 (15 years ago)

Author:

mido

Message:

a few moves of code from h to cpp, however, only part of the whole library is done

Location:

library

Files:

• bdm/base/bdmbase.cpp (modified) (5 diffs)
• bdm/base/bdmbase.h (modified) (8 diffs)
• bdm/base/datasources.cpp (modified) (1 diff)
• bdm/base/datasources.h (modified) (2 diffs)
• bdm/design/ctrlbase.cpp (modified) (2 diffs)
• bdm/design/ctrlbase.h (modified) (2 diffs)
• bdm/estim/arx.cpp (modified) (2 diffs)
• bdm/estim/arx.h (modified) (1 diff)
• bdm/math/chmat.cpp (modified) (1 diff)
• bdm/math/chmat.h (modified) (1 diff)
• tests/testsuite/randun_test.cpp (modified) (1 diff)

library/bdm/base/bdmbase.cpp

@@ -251 +251 @@
 }
 
+RV RV::expand_delayes() const {
+        RV rvt = this->remove_time(); //rv at t=0
+        RV tmp = rvt;
+        int td = mint();
+        for ( int i = -1; i >= td; i-- ) {
+                rvt.t_plus ( -1 );
+                tmp.add ( rvt ); //shift u1
+        }
+        return tmp;
+}
+
 str RV::tostr() const {
         ivec idlist ( dsize );
@@ -323 +334 @@
 }
 
+std::string RV::scalarname ( int scalat ) const {
+        bdm_assert ( scalat < dsize, "Wrong input index" );
+        int id = 0;
+        int scalid = 0;
+        while ( scalid + SIZES ( ids ( id ) ) <= scalat )  {
+                scalid += SIZES ( ids ( id ) );
+                id++;
+        };
+        //now id is the id of variable of interest
+        if ( size ( id ) == 1 )
+                return  NAMES ( ids ( id ) );
+        else
+                return  NAMES ( ids ( id ) ) + "_" + num2str ( scalat - scalid );
+
+}
+
 ivec RV::findself ( const RV &rv2 ) const {
         int i, j;
@@ -398 +425 @@
 }
 
+int logger::add_vector ( const RV &rv, string prefix ) {
+        int id;
+        if ( rv._dsize() > 0 ) {
+                id = entries.length();
+                names = concat ( names, prefix ); // diff
+                entries.set_length ( id + 1, true );
+                entries ( id ) = rv;
+        } else {
+                id = -1;
+        }
+        return id; // identifier of the last entry
+}
+
+int logger::add_setting ( const string &prefix ) {
+        Setting &root = setting_conf.getRoot();
+        int id = root.getLength(); //root must be group!!
+        if ( prefix.length() > 0 ) {
+                settings.set_length ( id + 1, true );
+                settings ( id ) = &root.add ( prefix, Setting::TypeList );
+        } else {
+                id = -1;
+        }
+        return id;
+}
+
+void epdf::log_register ( logger &L, const string &prefix ) {
+        RV r;
+        if ( isnamed() ) {
+                r = _rv();
+        } else {
+                r = RV ( "", dimension() );
+        };
+        root::log_register ( L, prefix );
+
+        // log full data
+        if ( log_level == 10 ) {
+                logrec->ids.set_size ( 1 );
+                logrec->ids ( 0 ) = logrec->L.add_setting ( prefix );
+        } else {
+                // log only
+                logrec->ids.set_size ( 3 );
+                if ( log_level > 0 ) {
+                        logrec->ids ( 0 ) = logrec->L.add_vector ( r, prefix + logrec->L.prefix_sep() + "mean" );
+                }
+                if ( log_level > 1 ) {
+                        logrec->ids ( 1 ) = logrec->L.add_vector ( r, prefix + logrec->L.prefix_sep() + "lb" );
+                        logrec->ids ( 2 ) = logrec->L.add_vector ( r, prefix + logrec->L.prefix_sep() + "ub" );
+                }
+        }
+}
+
+void epdf::log_write() const {
+        if ( log_level == 10 ) {
+                to_setting ( logrec->L.log_to_setting ( logrec->ids ( 0 ) ) );
+        } else {
+                if ( log_level > 0 ) {
+                        logrec->L.log_vector ( logrec->ids ( 0 ), mean() );
+                }
+                if ( log_level > 1 ) {
+                        vec lb;
+                        vec ub;
+                        qbounds ( lb, ub );
+                        logrec->L.log_vector ( logrec->ids ( 1 ), lb );
+                        logrec->L.log_vector ( logrec->ids ( 2 ), ub );
+                }
+        }
+}
+
+void datalink_buffered::set_connection ( const RV &rv, const RV &rv_up ) {
+        // create link between up and down
+        datalink_part::set_connection ( rv, rv_up.remove_time() ); // only non-delayed version
+
+        // create rvs of history
+        // we can store only what we get in rv_up - everything else is removed
+        ivec valid_ids = rv.findself_ids ( rv_up );
+        RV rv_hist = rv.subselect ( find ( valid_ids >= 0 ) ); // select only rvs that are in rv_up
+        RV rv_hist0 = rv_hist.remove_time(); // these RVs will form history at time =0
+        // now we need to know what is needed from Up
+        rv_hist = rv_hist.expand_delayes(); // full regressor - including time 0
+        Hrv = rv_hist.subt ( rv_hist0 );   // remove time 0
+        history = zeros ( Hrv._dsize() );
+
+        // decide if we need to copy val to history
+        if ( Hrv._dsize() > 0 ) {
+                v2h_up = rv_hist0.dataind ( rv_up ); // indeces of elements of rv_up to be copied
+        } // else v2h_up is empty
+
+        Hrv.dataind ( rv, h2v_hist, h2v_down );
+
+        downsize = v2v_down.length() + h2v_down.length();
+        upsize = v2v_up.length();
+}
+
+void datalink_buffered::set_history ( const RV& rv1, const vec &hist0 ) {
+        bdm_assert ( rv1._dsize() == hist0.length(), "hist is not compatible with given rv1" );
+        ivec ind_H;
+        ivec ind_h0;
+        Hrv.dataind ( rv1, ind_H, ind_h0 ); // find indeces of rv in
+        set_subvector ( history, ind_H, hist0 ( ind_h0 ) ); // copy given hist to appropriate places
+}
+
+void DS::log_register ( logger &L,  const string &prefix ) {
+        bdm_assert ( ytsize == Yrv._dsize(), "invalid DS: ytsize (" + num2str ( ytsize ) + ") different from Drv " + num2str ( Yrv._dsize() ) );
+        bdm_assert ( utsize == Urv._dsize(), "invalid DS: utsize (" + num2str ( utsize ) + ") different from Urv " + num2str ( Urv._dsize() ) );
+
+        root::log_register ( L, prefix );
+        //we know that
+        if ( log_level > 0 ) {
+                logrec->ids.set_size ( 2 );
+                logrec->ids ( 0 ) = logrec->L.add_vector ( Yrv, prefix );
+                logrec->ids ( 1 ) = logrec->L.add_vector ( Urv, prefix );
+        }
+}
+
+void DS::log_write ( ) const {
+        if ( log_level > 0 ) {
+                vec tmp ( Yrv._dsize() + Urv._dsize() );
+                getdata ( tmp );
+                // d is first in getdata
+                logrec->L.log_vector ( logrec->ids ( 0 ), tmp.left ( Yrv._dsize() ) );
+                // u follows after d in getdata
+                logrec->L.log_vector ( logrec->ids ( 1 ), tmp.mid ( Yrv._dsize(), Urv._dsize() ) );
+        }
+}
+
+void BM::set_options ( const string &opt ) {
+        if ( opt.find ( "logfull" ) != string::npos ) {
+                const_cast<epdf&> ( posterior() ).set_log_level ( 10 ) ;
+        } else {
+                if ( opt.find ( "logbounds" ) != string::npos ) {
+                        const_cast<epdf&> ( posterior() ).set_log_level ( 2 ) ;
+                } else {
+                        const_cast<epdf&> ( posterior() ).set_log_level ( 1 ) ;
+                }
+                if ( opt.find ( "logll" ) != string::npos ) {
+                        log_level = 1;
+                }
+        }
+}
+
+void BM::log_register ( logger &L, const string &prefix ) {
+        root::log_register ( L, prefix );
+
+        const_cast<epdf&> ( posterior() ).log_register ( L, prefix + L.prefix_sep() + "apost" );
+
+        if ( ( log_level ) > 0 ) {
+                logrec->ids.set_size ( 1 );
+                logrec->ids ( 0 ) = L.add_vector ( RV ( "", 1 ), prefix + L.prefix_sep() + "ll" );
+        }
+}
+
+void BM::log_write ( ) const {
+        posterior().log_write();
+        if ( log_level > 0 ) {
+                logrec->L.logit ( logrec->ids ( 0 ), ll );
+        }
+}
+
 void BM::bayes_batch ( const mat &Data, const vec &cond ) {
         for ( int t = 0; t < Data.cols(); t++ ) {
@@ -403 +588 @@
         }
 }
+
 void BM::bayes_batch ( const mat &Data, const mat &Cond ) {
         for ( int t = 0; t < Data.cols(); t++ ) {
@@ -408 +594 @@
         }
 }
-}
+
+}

library/bdm/base/bdmbase.h

@@ -186 +186 @@
         }
         //! returns name of a scalar at position scalat, i.e. it can be in the middle of vector name, in that case it adds "_%d" to it
-        std::string scalarname ( int scalat ) const {
-                bdm_assert ( scalat < dsize, "Wrong input index" );
-                int id = 0;
-                int scalid = 0;
-                while ( scalid + SIZES ( ids ( id ) ) <= scalat )  {
-                        scalid += SIZES ( ids ( id ) );
-                        id++;
-                };
-                //now id is the id of variable of interest
-                if ( size ( id ) == 1 )
-                        return  NAMES ( ids ( id ) );
-                else
-                        return  NAMES ( ids ( id ) ) + "_" + num2str ( scalat - scalid );
-
-        }
+        std::string scalarname ( int scalat ) const;
+
         void set_time ( int at, int time0 ) {
                 times ( at ) = time0;
@@ -246 +233 @@
         }
         //! return rvs with expanded delayes and sorted in the order of: \f$ [ rv_{0}, rv_{-1},\ldots  rv_{max_delay}]\f$
-        RV expand_delayes() const {
-                RV rvt = this->remove_time(); //rv at t=0
-                RV tmp = rvt;
-                int td = mint();
-                for ( int i = -1; i >= td; i-- ) {
-                        rvt.t_plus ( -1 );
-                        tmp.add ( rvt ); //shift u1
-                }
-                return tmp;
-        }
+        RV expand_delayes() const;
         //!@}
 
@@ -336 +314 @@
         //! returns an identifier which will be later needed for calling the \c logit() function
         //! For empty RV it returns -1, this entry will be ignored by \c logit().
-        virtual int add_vector ( const RV &rv, string prefix = "" ) {
-                int id;
-                if ( rv._dsize() > 0 ) {
-                        id = entries.length();
-                        names = concat ( names, prefix ); // diff
-                        entries.set_length ( id + 1, true );
-                        entries ( id ) = rv;
-                } else {
-                        id = -1;
-                }
-                return id; // identifier of the last entry
-        }
-        virtual int add_setting ( const string &prefix ) {
-                Setting &root = setting_conf.getRoot();
-                int id = root.getLength(); //root must be group!!
-                if ( prefix.length() > 0 ) {
-                        settings.set_length ( id + 1, true );
-                        settings ( id ) = &root.add ( prefix, Setting::TypeList );
-                } else {
-                        id = -1;
-                }
-                return id;
-        }
+        virtual int add_vector ( const RV &rv, string prefix = "" );
+
+        virtual int add_setting ( const string &prefix );
 
         //! log this vector
@@ -622 +580 @@
         //!  #1 mean
         //!  #2 mean + lower & upper bound
-        void log_register ( logger &L, const string &prefix ) {
-                RV r;
-                if ( isnamed() ) {
-                        r = _rv();
-                } else {
-                        r = RV ( "", dimension() );
-                };
-                root::log_register ( L, prefix );
-
-                // log full data
-                if ( log_level == 10 ) {
-                        logrec->ids.set_size ( 1 );
-                        logrec->ids ( 0 ) = logrec->L.add_setting ( prefix );
-                } else {
-                        // log only
-                        logrec->ids.set_size ( 3 );
-                        if ( log_level > 0 ) {
-                                logrec->ids ( 0 ) = logrec->L.add_vector ( r, prefix + logrec->L.prefix_sep() + "mean" );
-                        }
-                        if ( log_level > 1 ) {
-                                logrec->ids ( 1 ) = logrec->L.add_vector ( r, prefix + logrec->L.prefix_sep() + "lb" );
-                                logrec->ids ( 2 ) = logrec->L.add_vector ( r, prefix + logrec->L.prefix_sep() + "ub" );
-                        }
-                }
-        }
-        void log_write() const {
-                if ( log_level == 10 ) {
-                        to_setting ( logrec->L.log_to_setting ( logrec->ids ( 0 ) ) );
-                } else {
-                        if ( log_level > 0 ) {
-                                logrec->L.log_vector ( logrec->ids ( 0 ), mean() );
-                        }
-                        if ( log_level > 1 ) {
-                                vec lb;
-                                vec ub;
-                                qbounds ( lb, ub );
-                                logrec->L.log_vector ( logrec->ids ( 1 ), lb );
-                                logrec->L.log_vector ( logrec->ids ( 2 ), ub );
-                        }
-                }
-        }
+        void log_register ( logger &L, const string &prefix );
+
+        void log_write() const;
         //!@}
 
@@ -868 +788 @@
         }
 
-        void set_connection ( const RV &rv, const RV &rv_up ) {
-                // create link between up and down
-                datalink_part::set_connection ( rv, rv_up.remove_time() ); // only non-delayed version
-
-                // create rvs of history
-                // we can store only what we get in rv_up - everything else is removed
-                ivec valid_ids = rv.findself_ids ( rv_up );
-                RV rv_hist = rv.subselect ( find ( valid_ids >= 0 ) ); // select only rvs that are in rv_up
-                RV rv_hist0 = rv_hist.remove_time(); // these RVs will form history at time =0
-                // now we need to know what is needed from Up
-                rv_hist = rv_hist.expand_delayes(); // full regressor - including time 0
-                Hrv = rv_hist.subt ( rv_hist0 );   // remove time 0
-                history = zeros ( Hrv._dsize() );
-
-                // decide if we need to copy val to history
-                if ( Hrv._dsize() > 0 ) {
-                        v2h_up = rv_hist0.dataind ( rv_up ); // indeces of elements of rv_up to be copied
-                } // else v2h_up is empty
-
-                Hrv.dataind ( rv, h2v_hist, h2v_down );
-
-                downsize = v2v_down.length() + h2v_down.length();
-                upsize = v2v_up.length();
-        }
+        void set_connection ( const RV &rv, const RV &rv_up );
+        
         //! set history of variable given by \c rv1 to values of \c hist.
-        void set_history ( const RV& rv1, const vec &hist0 ) {
-                bdm_assert ( rv1._dsize() == hist0.length(), "hist is not compatible with given rv1" );
-                ivec ind_H;
-                ivec ind_h0;
-                Hrv.dataind ( rv1, ind_H, ind_h0 ); // find indeces of rv in
-                set_subvector ( history, ind_H, hist0 ( ind_h0 ) ); // copy given hist to appropriate places
-        }
+        void set_history ( const RV& rv1, const vec &hist0 );
 };
 
@@ -1060 +952 @@
 
         //! Register DS for logging into logger L
-        virtual void log_register ( logger &L,  const string &prefix ) {
-                bdm_assert ( ytsize == Yrv._dsize(), "invalid DS: ytsize (" + num2str ( ytsize ) + ") different from Drv " + num2str ( Yrv._dsize() ) );
-                bdm_assert ( utsize == Urv._dsize(), "invalid DS: utsize (" + num2str ( utsize ) + ") different from Urv " + num2str ( Urv._dsize() ) );
-
-                root::log_register ( L, prefix );
-                //we know that
-                if ( log_level > 0 ) {
-                        logrec->ids.set_size ( 2 );
-                        logrec->ids ( 0 ) = logrec->L.add_vector ( Yrv, prefix );
-                        logrec->ids ( 1 ) = logrec->L.add_vector ( Urv, prefix );
-                }
-        }
+        virtual void log_register ( logger &L,  const string &prefix );
         //! Register DS for logging into logger L
-        virtual void log_write ( ) const {
-                if ( log_level > 0 ) {
-                        vec tmp ( Yrv._dsize() + Urv._dsize() );
-                        getdata ( tmp );
-                        // d is first in getdata
-                        logrec->L.log_vector ( logrec->ids ( 0 ), tmp.left ( Yrv._dsize() ) );
-                        // u follows after d in getdata
-                        logrec->L.log_vector ( logrec->ids ( 1 ), tmp.mid ( Yrv._dsize(), Urv._dsize() ) );
-                }
-        }
+        virtual void log_write ( ) const;
         //!access function
         virtual const RV& _drv() const {
@@ -1245 +1117 @@
 
         //! Set boolean options from a string, recognized are: "logbounds,logll"
-        virtual void set_options ( const string &opt ) {
-                if ( opt.find ( "logfull" ) != string::npos ) {
-                        const_cast<epdf&> ( posterior() ).set_log_level ( 10 ) ;
-                } else {
-                        if ( opt.find ( "logbounds" ) != string::npos ) {
-                                const_cast<epdf&> ( posterior() ).set_log_level ( 2 ) ;
-                        } else {
-                                const_cast<epdf&> ( posterior() ).set_log_level ( 1 ) ;
-                        }
-                        if ( opt.find ( "logll" ) != string::npos ) {
-                                log_level = 1;
-                        }
-                }
-        }
+        virtual void set_options ( const string &opt );
 
         //! Add all logged variables to a logger
@@ -1264 +1123 @@
         //!  * y = 0/1 log-likelihood is to be logged
         //!  * x = level of the posterior (typically 0/1/2 for nothing/mean/bounds)
-        virtual void log_register ( logger &L, const string &prefix = "" ) {
-                root::log_register ( L, prefix );
-
-                const_cast<epdf&> ( posterior() ).log_register ( L, prefix + L.prefix_sep() + "apost" );
-
-                if ( ( log_level ) > 0 ) {
-                        logrec->ids.set_size ( 1 );
-                        logrec->ids ( 0 ) = L.add_vector ( RV ( "", 1 ), prefix + L.prefix_sep() + "ll" );
-                }
-        }
+        virtual void log_register ( logger &L, const string &prefix = "" );
+
         //! Save results to the given logger, details of what is stored is configured by \c LIDs and \c options
-        virtual void log_write ( ) const {
-                posterior().log_write();
-                if ( log_level > 0 ) {
-                        logrec->L.logit ( logrec->ids ( 0 ), ll );
-                }
-        }
+        virtual void log_write ( ) const;
+
         //!@}
         void from_setting ( const Setting &set ) {

library/bdm/base/datasources.cpp

@@ -40 +40 @@
         time = 0;
         Data = Dat;
+}
+
+void PdfDS::step() {
+        yt2rgr.store_data ( yt ); // y is now history
+        ut2rgr.filldown ( ut, rgr );
+        yt2rgr.filldown ( yt, rgr );
+        yt = ipdf->samplecond ( rgr );
+        ut2rgr.store_data ( ut ); //u is now history
+}
+
+void PdfDS::getdata ( vec &dt_out ) const {
+        bdm_assert_debug ( dt_out.length() >= utsize + ytsize, "Short output vector" );
+        dt_out.set_subvector ( 0, yt );
+        dt_out.set_subvector ( ytsize, ut );
+}
+
+void StateDS::step() {
+        vec imc ( IM->dimensionc() );
+        imc.set_subvector ( 0, xt );
+        u2imc.filldown ( ut, imc );
+        xt = IM->samplecond ( imc );
+
+        vec omc ( OM->dimensionc() );
+        omc.set_subvector ( 0, xt );
+        u2omc.filldown ( ut, omc );
+        vec yt;
+        yt = OM->samplecond ( omc );
+        //fill all data
+        dt.set_subvector ( 0, yt );
+        dt.set_subvector ( yt.length(), xt );
+        dt.set_subvector ( ytsize, ut );
 }
 

library/bdm/base/datasources.h

@@ -171 +171 @@
 
 public:
-        void step() {
-                yt2rgr.store_data ( yt ); // y is now history
-                ut2rgr.filldown ( ut, rgr );
-                yt2rgr.filldown ( yt, rgr );
-                yt = ipdf->samplecond ( rgr );
-                ut2rgr.store_data ( ut ); //u is now history
-        }
-        void getdata ( vec &dt_out ) const {
-                bdm_assert_debug ( dt_out.length() >= utsize + ytsize, "Short output vector" );
-                dt_out.set_subvector ( 0, yt );
-                dt_out.set_subvector ( ytsize, ut );
-        }
+        void step();
+
+        void getdata ( vec &dt_out ) const;
+
         void write ( const vec &ut0 ) {
                 ut = ut0;
@@ -350 +342 @@
         }
 
-        virtual void step() {
-                vec imc ( IM->dimensionc() );
-                imc.set_subvector ( 0, xt );
-                u2imc.filldown ( ut, imc );
-                xt = IM->samplecond ( imc );
-
-                vec omc ( OM->dimensionc() );
-                omc.set_subvector ( 0, xt );
-                u2omc.filldown ( ut, omc );
-                vec yt;
-                yt = OM->samplecond ( omc );
-                //fill all data
-                dt.set_subvector ( 0, yt );
-                dt.set_subvector ( yt.length(), xt );
-                dt.set_subvector ( ytsize, ut );
-        }
+        virtual void step();
 
         //! set parameters

library/bdm/design/ctrlbase.cpp

@@ -7 +7 @@
 }
 
+void LQG::update_system() {
+        pr.set_submatrix ( 0, 0, S->_B() );
+        pr.set_submatrix ( 0, dimu, S->_A() );
+
+        //penalization
+        qux.set_submatrix ( 0, 0, Qu._Ch() );
+        qux.set_submatrix ( 0, dimx + dimu + dimy, Qu._Ch() );
+
+        qyx.set_submatrix ( 0, 0, S->_C() );
+        qyx.set_submatrix ( 0, dimx, -eye ( dimy ) );
+
+        // parts of QR
+        hqy = Qy.to_mat() * qyx * pr;
+
+        // pre_qr
+        pre_qr = concat_vertical ( s * pr, concat_vertical ( hqy, qux ) );
+}
+        
 void LQG::set_control_parameters ( const mat &Qy0, const mat &Qu0, const vec &y_req0, int horizon0 ) {
         Qy = Qy0;
@@ -38 +56 @@
 }
 
+void LQG::ricatti_step() {
+        pre_qr.set_submatrix ( 0, 0, s*pr );
+        pre_qr.set_submatrix ( dimx + dimu + dimy, dimu + dimx, -Qy.to_mat() *y_req );
+        if ( !qr ( pre_qr, post_qr ) ) {
+                bdm_warning ( "QR in LQG unstable" );
+        }
+        triu ( post_qr );
+        // hn(m+1:2*m+n+r,m+1:2*m+n+r);
+        s = post_qr.get ( dimu, 2 * dimu + dimx + dimy - 1, dimu, 2 * dimu + dimx + dimy - 1 );
+};
+
+void LQG::redesign() {
+        for ( td = horizon; td > 0; td-- ) {
+                update_system();
+                ricatti_step();
+        }
+        /*                      ws=hn(1:m,m+1:2*m+n+r);
+                                wsd=hn(1:m,1:m);
+                                Lklq=-inv(wsd)*ws;*/
+        L = -inv ( post_qr.get ( 0, dimu - 1, 0, dimu - 1 ) ) * post_qr.get ( 0, dimu - 1, dimu, 2 * dimu + dimx + dimy - 1 );
+}
+
 
 }

library/bdm/design/ctrlbase.h

@@ -106 +106 @@
         void set_system ( shared_ptr<StateSpace<chmat> > S0 );
         //! update internal whan system has changed
-        void update_system() {
-                pr.set_submatrix ( 0, 0, S->_B() );
-                pr.set_submatrix ( 0, dimu, S->_A() );
-
-                //penalization
-                qux.set_submatrix ( 0, 0, Qu._Ch() );
-                qux.set_submatrix ( 0, dimx + dimu + dimy, Qu._Ch() );
-
-                qyx.set_submatrix ( 0, 0, S->_C() );
-                qyx.set_submatrix ( 0, dimx, -eye ( dimy ) );
-
-                // parts of QR
-                hqy = Qy.to_mat() * qyx * pr;
-
-                // pre_qr
-                pre_qr = concat_vertical ( s * pr, concat_vertical ( hqy, qux ) );
-        }
+        void update_system();
         //! set penalization matrices and control horizon
         void set_control_parameters ( const mat &Qy0, const mat &Qu0, const vec &y_req0, int horizon0 );
@@ -137 +121 @@
         //! function for future use which is called at each time td; Should call initialize()!
         //! redesign one step of the
-        void ricatti_step() {
-                pre_qr.set_submatrix ( 0, 0, s*pr );
-                pre_qr.set_submatrix ( dimx + dimu + dimy, dimu + dimx, -Qy.to_mat() *y_req );
-                if ( !qr ( pre_qr, post_qr ) ) {
-                        bdm_warning ( "QR in LQG unstable" );
-                }
-                triu ( post_qr );
-                // hn(m+1:2*m+n+r,m+1:2*m+n+r);
-                s = post_qr.get ( dimu, 2 * dimu + dimx + dimy - 1, dimu, 2 * dimu + dimx + dimy - 1 );
-        };
-        void redesign() {
-                for ( td = horizon; td > 0; td-- ) {
-                        update_system();
-                        ricatti_step();
-                }
-                /*                      ws=hn(1:m,m+1:2*m+n+r);
-                                        wsd=hn(1:m,1:m);
-                                        Lklq=-inv(wsd)*ws;*/
-                L = -inv ( post_qr.get ( 0, dimu - 1, 0, dimu - 1 ) ) * post_qr.get ( 0, dimu - 1, dimu, 2 * dimu + dimx + dimy - 1 );
-        }
+        void ricatti_step();
+        
+        void redesign();
+
         //! compute control action
         vec ctrlaction ( const vec &state, const vec &ukm ) const {

library/bdm/estim/arx.cpp

@@ -68 +68 @@
 }
 
+void ARX::flatten ( const BMEF* B ) {
+        const ARX* A = dynamic_cast<const ARX*> ( B );
+        // nu should be equal to B.nu
+        est.pow ( A->posterior()._nu() / posterior()._nu() );
+        if ( evalll ) {
+                last_lognc = est.lognc();
+        }
+}
+
 ARX* ARX::_copy_ ( ) const {
         ARX* Tmp = new ARX ( *this );
@@ -98 +107 @@
 }
 
+enorm<ldmat>* ARX::epredictor() const {
+        bdm_assert_debug ( dimy == posterior()._V().rows() - 1, "Regressor is not only 1" );
+        return epredictor ( vec_1 ( 1.0 ) );
+}
 
 mlstudent* ARX::predictor_student ( ) const {

library/bdm/estim/arx.h

@@ -82 +82 @@
         };
         double logpred ( const vec &yt ) const;
-        void flatten ( const BMEF* B ) {
-                const ARX* A = dynamic_cast<const ARX*> ( B );
-                // nu should be equal to B.nu
-                est.pow ( A->posterior()._nu() / posterior()._nu() );
-                if ( evalll ) {
-                        last_lognc = est.lognc();
-                }
-        }
+        void flatten ( const BMEF* B );
         //! Conditioned version of the predictor
         enorm<ldmat>* epredictor ( const vec &rgr ) const;
         //! Predictor for empty regressor
-        enorm<ldmat>* epredictor() const {
-                bdm_assert_debug ( dimy == posterior()._V().rows() - 1, "Regressor is not only 1" );
-                return epredictor ( vec_1 ( 1.0 ) );
-        }
+        enorm<ldmat>* epredictor() const;
         //! conditional version of the predictor
         template<class sq_T>

library/bdm/math/chmat.cpp

@@ -5 +5 @@
 
 namespace bdm {
+
+void chmat::add ( const chmat &A2, double w ) {
+        bdm_assert_debug ( dim == A2.dim, "Matrices of unequal dimension" );
+        mat pre = concat_vertical ( Ch, sqrt ( w ) * A2.Ch );
+        mat post = zeros ( pre.rows(), pre.cols() );
+        if ( !qr ( pre, post ) ) {
+                bdm_warning ( "Unstable QR in chmat add" );
+        }
+        Ch = post ( 0, dim - 1, 0, dim - 1 );
+};
 
 void chmat::opupdt ( const vec &v, double w ) {

library/bdm/math/chmat.h

@@ -44 +44 @@
         void clear();
         //! add another chmat \c A2 with weight \c w.
-        void add ( const chmat &A2, double w = 1.0 ) {
-                bdm_assert_debug ( dim == A2.dim, "Matrices of unequal dimension" );
-                mat pre = concat_vertical ( Ch, sqrt ( w ) * A2.Ch );
-                mat post = zeros ( pre.rows(), pre.cols() );
-                if ( !qr ( pre, post ) ) {
-                        bdm_warning ( "Unstable QR in chmat add" );
-                }
-                Ch = post ( 0, dim - 1, 0, dim - 1 );
-        };
+        void add ( const chmat &A2, double w = 1.0 );
         //!Inversion in the same form, i.e. cholesky
         void inv ( chmat &Inv ) const   {

library/tests/testsuite/randun_test.cpp

@@ -5 +5 @@
 using namespace itpp;
 
-TEST ( test_randun ) {
+TEST ( randun_test ) {
         // matlab output obtained by
         // >>clear all