Changeset 598

Show
Ignore:
Timestamp:
09/03/09 00:27:23 (15 years ago)
Author:
smidl
Message:

new buffered datalink ( #32 ) and new datasources - all with minor trivial tests

Location:
library
Files:
3 added
7 modified

Legend:

Unmodified
Added
Removed

  • library/bdm/base/bdmbase.cpp

    r565 r598  
    191191void datalink::set_connection (const RV &rv, const RV &rv_up) { 
    192192        downsize = rv._dsize(); 
    193         upsize = rv_up._dsize(); 
    194         v2v_up = rv.dataind (rv_up); 
    195  
    196         bdm_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up"); 
     193          upsize = rv_up._dsize(); 
     194          v2v_up = rv.dataind (rv_up); 
     195          bdm_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up"); 
    197196} 
    198197 
    … …  
    201200        upsize = us; 
    202201        v2v_up = upind; 
    203  
    204         bdm_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up"); 
     202    bdm_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up"); 
     203} 
     204 
     205void datalink_part::set_connection (const RV &rv, const RV &rv_up) { 
     206         rv.dataind (rv_up, v2v_down, v2v_up); 
     207         downsize = v2v_down.length(); 
     208         upsize = v2v_up.length(); 
    205209} 
    206210 
    … …  
    321325} 
    322326 
     327ivec RV::findself_ids ( const RV &rv2 ) const { 
     328        int i, j; 
     329        ivec tmp = -ones_i ( len ); 
     330        for ( i = 0; i < len; i++ ) { 
     331                for ( j = 0; j < rv2.length(); j++ ) { 
     332                        if ( ( ids ( i ) == rv2.ids ( j ) ) ) { 
     333                                tmp ( i ) = j; 
     334                                break; 
     335                        } 
     336                } 
     337        } 
     338        return tmp; 
     339} 
     340 
    323341void RV::from_setting ( const Setting &set ) { 
    324342        Array<string> A; 

  • library/bdm/base/bdmbase.h

    r596 r598  
    9999 
    100100        private: 
     101          enum enum_dummy {dummy}; 
    101102                //! auxiliary function used in constructor 
    102103                void init (const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times); 
    103104                int init (const string &name, int size); 
     105                //! Private constructor from IDs, potentially dangerous since all ids must be valid! 
     106                //! dummy is there to prevent confusion with RV(" string"); 
     107                explicit RV (const ivec &ids0, enum_dummy dum): dsize(0), len(ids0.length()), ids(ids0), times(zeros_i(ids0.length())) { 
     108                  dsize = countsize(); 
     109                } 
    104110        public: 
    105111                //! \name Constructors 
    … …  
    171177                //! Find indices of self in another rv, \return ivec of the same size as self. 
    172178                ivec findself (const RV &rv2) const; 
     179                //! Find indices of self in another rv, ignore time, \return ivec of the same size as self. 
     180                ivec findself_ids (const RV &rv2) const; 
    173181                //! Compare if \c rv2 is identical to this \c RV 
    174182                bool equal (const RV &rv2) const; 
    … …  
    192200                //!@} 
    193201 
     202                //! @{ \name Auxiliary functions 
     203                //! returns rvs with time set to 0 and removed duplicates 
     204                RV remove_time() const {  return RV(unique(ids), dummy);        } 
     205                //! return rvs with expanded delayes and sorted in the order of: \f$ [ rv_{0}, rv_{-1}, rv_{ 
     206                RV expand_delayes() const { 
     207                  RV rvt=this->remove_time(); //rv at t=0 
     208                  RV tmp = rvt; 
     209                  int td = mint(); 
     210                  for ( int i = -1; i >= td; i-- ) { 
     211                        rvt.t ( -1 ); 
     212                        tmp.add ( rvt ); //shift u1 
     213                  } 
     214                  return tmp; 
     215                } 
     216                //!@} 
     217                 
    194218                //!\name Relation to vectors 
    195219                //!@{ 
    … …  
    200224                //! Then, data can be copied via: data_of_this = cdata(ind); 
    201225                ivec dataind (const RV &crv) const; 
     226                //! same as dataind but this time crv should not be complete supperset of rv. 
     227                ivec dataind_part (const RV &crv) const; 
    202228                //! generate mutual indices when copying data between self and crv. 
    203229                //! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i) 
    … …  
    461487                //! @{ 
    462488 
    463                 RV _rv() const { 
     489                const RV& _rv() const { 
    464490                        return ep->_rv(); 
    465491                } 
    466                 RV _rvc() { 
     492                const RV& _rvc() { 
    467493                        return rvc; 
    468494                } 
    … …  
    570596                //! val-to-val link, indices of the upper val 
    571597                ivec v2v_up; 
    572  
     598                 
    573599        public: 
    574600                //! Constructor 
    … …  
    588614                //! Get val for myself from val of "Up" 
    589615                vec pushdown (const vec &val_up) { 
     616                  vec tmp(downsize); 
     617                  filldown(val_up, tmp); 
     618                        return tmp; 
     619                } 
     620                //! Get val for vector val_down from val of "Up" 
     621                void filldown (const vec &val_up, vec &val_down) { 
    590622                        bdm_assert_debug (upsize == val_up.length(), "Wrong val_up"); 
    591                         return get_vec (val_up, v2v_up); 
     623                        val_down=val_up(v2v_up); 
    592624                } 
    593625 
    … …  
    598630                        set_subvector (val_up, v2v_up, val); 
    599631                } 
     632                //! access functions 
     633                int _upsize(){return upsize;} 
     634                //! access functions 
     635                int _downsize(){return downsize;} 
     636}; 
     637 
     638/*! Extension of datalink to fill only part of Down 
     639*/ 
     640class datalink_part : public datalink{ 
     641  protected: 
     642        //! indeces of values in vector downsize 
     643        ivec v2v_down; 
     644  public: 
     645        void set_connection(const RV &rv, const RV &rv_up); 
     646                //! Get val for vector val_down from val of "Up" 
     647                void filldown (const vec &val_up, vec &val_down) { 
     648                        set_subvector(val_down, v2v_down, val_up(v2v_up)); 
     649                } 
     650}; 
     651 
     652/*! \brief Datalink that buffers delayed values - do not forget to call step() 
     653 
     654Up is current data, Down is their subset with possibly delayed values 
     655*/ 
     656class datalink_buffered: public datalink_part{ 
     657  protected: 
     658        //! History, ordered as \f$[Up_{t-1},Up_{t-2}, \ldots]\f$ 
     659        vec history; 
     660        //! rv of the history 
     661        RV Hrv; 
     662        //! h2v : indeces in down  
     663        ivec h2v_down; 
     664        //! h2v : indeces in history 
     665        ivec h2v_hist; 
     666  public: 
     667         
     668    datalink_buffered():datalink_part(),history(0), h2v_down(0), h2v_hist(0){}; 
     669        //! push current data to history 
     670        void step(const vec &val_up){ 
     671                if (Hrv._dsize()>0){ 
     672                  history.shift_right ( 0, Hrv._dsize() ); 
     673                  history.set_subvector ( 0, val_up ); // write U after Drv 
     674                } 
     675        } 
     676                //! Get val for myself from val of "Up" 
     677                vec pushdown (const vec &val_up) { 
     678                  vec tmp(downsize); 
     679                  filldown(val_up, tmp); 
     680                        return tmp; 
     681                } 
     682         
     683        void filldown(const vec &val_up, vec &val_down){ 
     684          bdm_assert_debug(val_down.length()>=downsize, "short val_down"); 
     685         
     686          set_subvector(val_down, v2v_down, val_up(v2v_up)); // copy direct values 
     687          set_subvector(val_down, h2v_down, history(h2v_hist)); // copy delayed values 
     688        } 
     689         
     690        void set_connection(const RV &rv, const RV &rv_up){ 
     691          // create link between up and down 
     692          datalink_part::set_connection(rv, rv_up); 
     693           
     694          // create rvs of history  
     695          // we can store only what we get in rv_up - everything else is removed 
     696          ivec valid_ids=rv.findself_ids(rv_up); 
     697          RV rv_hist = rv.subselect( find(valid_ids>=0) ); // select only rvs that are in rv_up 
     698          rv_hist = rv_hist.expand_delayes(); // full regressor - including time 0 
     699          Hrv=rv_hist.subt(rv.remove_time());   // remove time 0 
     700          history = zeros(Hrv._dsize()); 
     701           
     702          Hrv.dataind(rv,h2v_hist,h2v_down); 
     703           
     704          downsize = v2v_down.length()+h2v_down.length(); 
     705          upsize = v2v_up.length(); 
     706        } 
     707}; 
     708 
     709//! buffered datalink from 2 vectors to 1 
     710class datalink_2to1_buffered{ 
     711  protected: 
     712        datalink_buffered dl1; 
     713        datalink_buffered dl2; 
     714  public: 
     715        //! set connection between RVs 
     716        void set_connection(const RV &rv, const RV &rv_up1, const RV &rv_up2){ 
     717          dl1.set_connection(rv,rv_up1); 
     718          dl2.set_connection(rv,rv_up2); 
     719        } 
     720        void filldown(const vec &val1, const vec &val2, vec &val_down){ 
     721          bdm_assert_debug(val_down.length()>=dl1._downsize()+dl2._downsize(), "short val_down"); 
     722          dl1.filldown(val1, val_down); 
     723          dl2.filldown(val2, val_down); 
     724        } 
     725        void step(const vec &dt, const vec &ut) {dl1.step(dt); dl2.step(ut);} 
    600726}; 
    601727 
    … …  
    805931                } 
    806932                //!access function 
    807                 virtual RV _drv() const { 
     933                virtual const RV& _drv() const { 
    808934        //              return concat (Drv, Urv);// why!!! 
    809935                        return Drv;// why!!! 

  • library/bdm/base/datasources.h

    r565 r598  
    2828*/ 
    2929class MemDS : public DS { 
    30 protected: 
    31         //! internal matrix of data 
    32         mat Data; 
    33         //! active column in the Data matrix 
    34         int time; 
    35         //!  vector of rows that are presented in Dt 
    36         ivec rowid; 
    37         //! vector of delays that are presented in Dt 
    38         ivec delays; 
    39  
    40 public: 
    41         void getdata ( vec &dt ); 
    42         void getdata ( vec &dt, const ivec &indeces ); 
    43         void set_rvs ( RV &drv, RV &urv ); 
    44  
    45         void write ( vec &ut ) { 
    46                 bdm_error ( "MemDS::write is not supported" ); 
    47         } 
    48  
    49         void write ( vec &ut, ivec &indices ) { 
    50                 bdm_error ( "MemDS::write is not supported" ); 
    51         } 
    52  
    53         void step(); 
    54         //!Default constructor 
    55         MemDS () {}; 
    56         MemDS ( mat &Dat, ivec &rowid0, ivec &delays0 ); 
    57 }; 
     30        protected: 
     31                //! internal matrix of data 
     32                mat Data; 
     33                //! active column in the Data matrix 
     34                int time; 
     35                //!  vector of rows that are presented in Dt 
     36                ivec rowid; 
     37                //! vector of delays that are presented in Dt 
     38                ivec delays; 
     39 
     40        public: 
     41                void getdata ( vec &dt ); 
     42                void getdata ( vec &dt, const ivec &indeces ); 
     43                void set_rvs ( RV &drv, RV &urv ); 
     44 
     45                void write ( vec &ut ) { 
     46                        bdm_error ( "MemDS::write is not supported" ); 
     47                } 
     48 
     49                void write ( vec &ut, ivec &indices ) { 
     50                        bdm_error ( "MemDS::write is not supported" ); 
     51                } 
     52 
     53                void step(); 
     54                //!Default constructor 
     55                MemDS () {}; 
     56                MemDS ( mat &Dat, ivec &rowid0, ivec &delays0 ); 
     57}; 
     58 
     59/*!  \brief Simulate data from a static pdf 
     60Trivial example of a data source, could be used for tests of some estimation algorithms. For example, simulating data from a mixture model and feeding them to mixture model estimators. 
     61*/ 
     62 
     63class EpdfDS: public DS { 
     64        protected: 
     65                //! internal pointer to epdf from which we samplecond 
     66                shared_ptr<epdf> iepdf; 
     67                //! internal storage of data sample 
     68                vec dt; 
     69        public: 
     70                void step() { 
     71                        dt=iepdf->sample(); 
     72                } 
     73                void getdata ( vec &dt_out ) { 
     74                        dt_out = dt; 
     75                } 
     76                void getdata ( vec &dt_out, const ivec &ids ) { 
     77                        dt_out = dt ( ids ); 
     78                } 
     79                const RV& _drv() { 
     80                        return iepdf->_rv(); 
     81                } 
     82 
     83                /*! 
     84                \code 
     85                class = "PdfDS"; 
     86                epdf = {class="epdf_offspring", ...}// list of points 
     87                \endcode 
     88 
     89                */ 
     90                void from_setting ( const Setting &set ) { 
     91                        iepdf=UI::build<epdf> ( set,"epdf",UI::compulsory ); 
     92                        dt = zeros ( iepdf->dimension() ); 
     93                } 
     94}; 
     95UIREGISTER ( EpdfDS ); 
     96 
     97/*!  \brief Simulate data from conditional density 
     98Still having only one density but allowing conditioning on either input or delayed values. 
     99*/ 
     100class MpdfDS :public DS { 
     101        protected: 
     102                //! internal pointer to epdf from which we samplecond 
     103                shared_ptr<mpdf> impdf; 
     104                //! internal storage of data sample 
     105                vec dt; 
     106                //! input vector 
     107                vec ut; 
     108                //! datalink between dt and ut and regressor 
     109                datalink_2to1_buffered rgrlink; 
     110                //! numeric values of regressor 
     111                vec rgr; 
     112                 
     113        public: 
     114                void step() { 
     115                        rgrlink.filldown ( dt,ut,rgr ); 
     116                        rgrlink.step(dt,ut);//whist history 
     117                        dt=impdf->samplecond ( rgr ); 
     118                } 
     119                void getdata ( vec &dt_out ) { 
     120                        dt_out = dt; 
     121                } 
     122                void getdata ( vec &dt_out, const ivec &ids ) { 
     123                        dt_out = dt ( ids ); 
     124                } 
     125                const RV& _drv() const { 
     126                        return impdf->_rv(); 
     127                } 
     128                void write(const vec &ut0){ut=ut0;} 
     129 
     130                /*! 
     131                \code 
     132                class = "MpdfDS"; 
     133                mpdf = {class="epdf_offspring", ...}// list of points 
     134                \endcode 
     135 
     136                */ 
     137                void from_setting ( const Setting &set ) { 
     138                        impdf=UI::build<mpdf> ( set,"mpdf",UI::compulsory ); 
     139                         
     140                        // get unique rvs form rvc 
     141                        RV rgrv0=impdf->_rvc().remove_time(); 
     142                        // input is what in not in _rv() 
     143                        RV urv=rgrv0.subt(impdf->_rv());  
     144                        set_drv(impdf->_rv(), urv); 
     145                        // connect input and output to rvc 
     146                        rgrlink.set_connection(impdf->_rvc(), Drv,Urv);  
     147 
     148                        dt = zeros ( impdf->dimension() ); 
     149                        rgr = zeros ( impdf->dimensionc() ); 
     150                        ut = zeros(urv._dsize()); 
     151                } 
     152}; 
     153UIREGISTER ( MpdfDS ); 
    58154 
    59155/*! Pseudovirtual class for reading data from files 
    … …  
    62158class FileDS: public MemDS { 
    63159 
    64 public: 
    65         void getdata ( vec &dt ) { 
    66                 dt = Data.get_col ( time ); 
    67         } 
    68  
    69         void getdata ( vec &dt, const ivec &indices ) { 
    70                 vec tmp = Data.get_col ( time ); 
    71                 dt = tmp ( indices ); 
    72         } 
    73  
    74         //! returns number of data in the file; 
    75         int ndat() { 
    76                 return Data.cols(); 
    77         } 
    78         //! no sense to log this type 
    79         void log_add ( logger &L ) {}; 
    80         //! no sense to log this type 
    81         void logit ( logger &L ) {}; 
     160        public: 
     161                void getdata ( vec &dt ) { 
     162                        dt = Data.get_col ( time ); 
     163                } 
     164 
     165                void getdata ( vec &dt, const ivec &indices ) { 
     166                        vec tmp = Data.get_col ( time ); 
     167                        dt = tmp ( indices ); 
     168                } 
     169 
     170                //! returns number of data in the file; 
     171                int ndat() { 
     172                        return Data.cols(); 
     173                } 
     174                //! no sense to log this type 
     175                void log_add ( logger &L ) {}; 
     176                //! no sense to log this type 
     177                void logit ( logger &L ) {}; 
    82178}; 
    83179 
    … …  
    90186class ITppFileDS: public FileDS { 
    91187 
    92 public: 
    93         ITppFileDS ( const string &fname, const string &varname ) : FileDS() { 
    94                 it_file it ( fname ); 
    95                 it << Name ( varname ); 
    96                 it >> Data; 
    97                 time = 0; 
    98                 //rowid and delays are ignored 
    99         }; 
    100  
    101         ITppFileDS () : FileDS() { 
    102         }; 
    103  
    104         void from_setting ( const Setting &set ); 
    105  
    106         // TODO dodelat void to_setting( Setting &set ) const; 
     188        public: 
     189                ITppFileDS ( const string &fname, const string &varname ) : FileDS() { 
     190                        it_file it ( fname ); 
     191                        it << Name ( varname ); 
     192                        it >> Data; 
     193                        time = 0; 
     194                        //rowid and delays are ignored 
     195                }; 
     196 
     197                ITppFileDS () : FileDS() { 
     198                }; 
     199 
     200                void from_setting ( const Setting &set ); 
     201 
     202                // TODO dodelat void to_setting( Setting &set ) const; 
    107203 
    108204}; 
    … …  
    120216class CsvFileDS: public FileDS { 
    121217 
    122 public: 
    123         //! Constructor - create DS from a CSV file. 
    124         CsvFileDS ( const string& fname, const string& orientation = "BY_COL" ); 
     218        public: 
     219                //! Constructor - create DS from a CSV file. 
     220                CsvFileDS ( const string& fname, const string& orientation = "BY_COL" ); 
    125221}; 
    126222 
    … …  
    132228*/ 
    133229class ArxDS : public DS { 
    134 protected: 
    135         //! Rv of the regressor 
    136         RV Rrv; 
    137         //! History, ordered as \f$[y_t, u_t, y_{t-1 }, u_{t-1}, \ldots]\f$ 
    138         vec H; 
    139         //! (future) input 
    140         vec U; 
    141         //! temporary variable for regressor 
    142         vec rgr; 
    143         //! data link: H -> rgr 
    144         datalink rgrlnk; 
    145         //! model of Y - linear Gaussian 
    146         mlnorm<chmat> model; 
    147         //! options 
    148         bool opt_L_theta; 
    149         //! loggers 
    150         int L_theta; 
    151         int L_R; 
    152         int dt_size; 
    153 public: 
    154         void getdata ( vec &dt ) { 
    155                 dt = H; 
    156         } 
    157  
    158         void getdata ( vec &dt, const ivec &indices ) { 
    159                 dt = H ( indices ); 
    160         } 
    161  
    162         void write ( vec &ut ) { 
    163                 U = ut; 
    164         } 
    165  
    166         void write ( vec &ut, const ivec &indices ) { 
    167                 bdm_assert_debug ( ut.length() == indices.length(), "ArxDS" ); 
    168                 set_subvector ( U, indices, ut ); 
    169         } 
    170  
    171         void step(); 
    172  
    173         //!Default constructor 
    174         ArxDS ( ) {}; 
    175         //! Set parameters of the internal model, H is maximum time delay 
    176         void set_parameters ( const mat &Th0, const vec mu0, const chmat &sqR0 ) { 
    177                 model.set_parameters ( Th0, mu0, sqR0 ); 
    178         }; 
    179         //! Set 
    180         void set_drv ( const RV &yrv, const RV &urv, const RV &rrv ) { 
    181                 Rrv = rrv; 
    182                 Urv = urv; 
    183                 dt_size = yrv._dsize() + urv._dsize(); 
    184  
    185                 RV drv = concat ( yrv, urv ); 
    186                 Drv = drv; 
    187                 int td = rrv.mint(); 
    188                 H.set_size ( drv._dsize() * ( -td + 1 ) ); 
    189                 U.set_size ( Urv._dsize() ); 
    190                 for ( int i = -1; i >= td; i-- ) { 
    191                         drv.t ( -1 ); 
    192                         Drv.add ( drv ); //shift u1 
    193                 } 
    194                 rgrlnk.set_connection ( rrv, Drv ); 
    195  
    196                 dtsize = Drv._dsize(); 
    197                 utsize = Urv._dsize(); 
    198         } 
    199         //! set options from a string 
    200         void set_options ( const string &s ) { 
    201                 opt_L_theta = ( s.find ( "L_theta" ) != string::npos ); 
    202         }; 
    203         virtual void log_add ( logger &L ) { 
    204                 //DS::log_add ( L ); too long!! 
    205                 L_dt = L.add ( Drv ( 0, dt_size ), "" ); 
    206                 L_ut = L.add ( Urv, "" ); 
    207  
    208                 mat &A = model._A(); 
    209                 mat R = model._R(); 
    210                 if ( opt_L_theta ) { 
    211                         L_theta = L.add ( RV ( "{th }", vec_1 ( A.rows() * A.cols() ) ), "t" ); 
    212                 } 
    213                 if ( opt_L_theta ) { 
    214                         L_R = L.add ( RV ( "{R }", vec_1 ( R.rows() * R.cols() ) ), "r" ); 
    215                 } 
    216         } 
    217         virtual void logit ( logger &L ) { 
    218                 //DS::logit ( L ); 
    219                 L.logit ( L_dt, H.left ( dt_size ) ); 
    220                 L.logit ( L_ut, U ); 
    221  
    222                 mat &A = model._A(); 
    223                 mat R = model._R(); 
    224                 if ( opt_L_theta ) { 
    225                         L.logit ( L_theta, vec ( A._data(), A.rows() *A.cols() ) ); 
    226                 }; 
    227                 if ( opt_L_theta ) { 
    228                         L.logit ( L_R, vec ( R._data(), R.rows() *R.rows() ) ); 
    229                 }; 
    230         } 
    231  
    232         // TODO dokumentace - aktualizovat 
    233         /*! UI for ArxDS using factorized description! 
    234  
    235         The ArxDS is constructed from a structure with fields: 
    236         \code 
    237         system = { 
    238                 type = "ArxDS"; 
    239                 // description of y variables 
    240                 y = {type="rv"; names=["y", "u"];}; 
    241                 // description of u variable 
    242                 u = {type="rv"; names=[];} 
    243                 // description of regressor 
    244                 rgr = {type="rv"; 
    245                         names = ["y","y","y","u"]; 
    246                         times = [-1, -2, -3, -1]; 
    247                 } 
    248  
    249                 // theta 
    250                 theta = [0.8, -0.3, 0.4, 1.0, 
    251                                 0.0, 0.0, 0.0, 0.0]; 
    252                 // offset (optional) 
    253                 offset = [0.0, 0.0]; 
    254                 //variance 
    255                 r = [0.1, 0.0, 
    256                         0.0, 1.0]; 
    257                 //options: L_theta = log value of theta, 
    258                 opt = "L_theta"; 
    259         }; 
    260         \endcode 
    261  
    262         Result is ARX data source offering with full history as Drv. 
    263         */ 
    264         void from_setting ( const Setting &set ); 
    265  
    266         // TODO dodelat void to_setting( Setting &set ) const; 
     230        protected: 
     231                //! Rv of the regressor 
     232                RV Rrv; 
     233                //! History, ordered as \f$[y_t, u_t, y_{t-1 }, u_{t-1}, \ldots]\f$ 
     234                vec H; 
     235                //! (future) input 
     236                vec U; 
     237                //! temporary variable for regressor 
     238                vec rgr; 
     239                //! data link: H -> rgr 
     240                datalink rgrlnk; 
     241                //! model of Y - linear Gaussian 
     242                mlnorm<chmat> model; 
     243                //! options 
     244                bool opt_L_theta; 
     245                //! loggers 
     246                int L_theta; 
     247                int L_R; 
     248                int dt_size; 
     249        public: 
     250                void getdata ( vec &dt ) { 
     251                        dt = H; 
     252                } 
     253 
     254                void getdata ( vec &dt, const ivec &indices ) { 
     255                        dt = H ( indices ); 
     256                } 
     257 
     258                void write ( vec &ut ) { 
     259                        U = ut; 
     260                } 
     261 
     262                void write ( vec &ut, const ivec &indices ) { 
     263                        bdm_assert_debug ( ut.length() == indices.length(), "ArxDS" ); 
     264                        set_subvector ( U, indices, ut ); 
     265                } 
     266 
     267                void step(); 
     268 
     269                //!Default constructor 
     270                ArxDS ( ) {}; 
     271                //! Set parameters of the internal model, H is maximum time delay 
     272                void set_parameters ( const mat &Th0, const vec mu0, const chmat &sqR0 ) { 
     273                        model.set_parameters ( Th0, mu0, sqR0 ); 
     274                }; 
     275                //! Set 
     276                void set_drv ( const RV &yrv, const RV &urv, const RV &rrv ) { 
     277                        Rrv = rrv; 
     278                        Urv = urv; 
     279                        dt_size = yrv._dsize() + urv._dsize(); 
     280 
     281                        RV drv = concat ( yrv, urv ); 
     282                        Drv = drv; 
     283                        int td = rrv.mint(); 
     284                        H.set_size ( drv._dsize() * ( -td + 1 ) ); 
     285                        U.set_size ( Urv._dsize() ); 
     286                        for ( int i = -1; i >= td; i-- ) { 
     287                                drv.t ( -1 ); 
     288                                Drv.add ( drv ); //shift u1 
     289                        } 
     290                        rgrlnk.set_connection ( rrv, Drv ); 
     291 
     292                        dtsize = Drv._dsize(); 
     293                        utsize = Urv._dsize(); 
     294                } 
     295                //! set options from a string 
     296                void set_options ( const string &s ) { 
     297                        opt_L_theta = ( s.find ( "L_theta" ) != string::npos ); 
     298                }; 
     299                virtual void log_add ( logger &L ) { 
     300                        //DS::log_add ( L ); too long!! 
     301                        L_dt = L.add ( Drv ( 0, dt_size ), "" ); 
     302                        L_ut = L.add ( Urv, "" ); 
     303 
     304                        mat &A = model._A(); 
     305                        mat R = model._R(); 
     306                        if ( opt_L_theta ) { 
     307                                L_theta = L.add ( RV ( "{th }", vec_1 ( A.rows() * A.cols() ) ), "t" ); 
     308                        } 
     309                        if ( opt_L_theta ) { 
     310                                L_R = L.add ( RV ( "{R }", vec_1 ( R.rows() * R.cols() ) ), "r" ); 
     311                        } 
     312                } 
     313                virtual void logit ( logger &L ) { 
     314                        //DS::logit ( L ); 
     315                        L.logit ( L_dt, H.left ( dt_size ) ); 
     316                        L.logit ( L_ut, U ); 
     317 
     318                        mat &A = model._A(); 
     319                        mat R = model._R(); 
     320                        if ( opt_L_theta ) { 
     321                                L.logit ( L_theta, vec ( A._data(), A.rows() *A.cols() ) ); 
     322                        }; 
     323                        if ( opt_L_theta ) { 
     324                                L.logit ( L_R, vec ( R._data(), R.rows() *R.rows() ) ); 
     325                        }; 
     326                } 
     327 
     328                // TODO dokumentace - aktualizovat 
     329                /*! UI for ArxDS using factorized description! 
     330 
     331                The ArxDS is constructed from a structure with fields: 
     332                \code 
     333                system = { 
     334                        type = "ArxDS"; 
     335                        // description of y variables 
     336                        y = {type="rv"; names=["y", "u"];}; 
     337                        // description of u variable 
     338                        u = {type="rv"; names=[];} 
     339                        // description of regressor 
     340                        rgr = {type="rv"; 
     341                                names = ["y","y","y","u"]; 
     342                                times = [-1, -2, -3, -1]; 
     343                        } 
     344 
     345                        // theta 
     346                        theta = [0.8, -0.3, 0.4, 1.0, 
     347                                        0.0, 0.0, 0.0, 0.0]; 
     348                        // offset (optional) 
     349                        offset = [0.0, 0.0]; 
     350                        //variance 
     351                        r = [0.1, 0.0, 
     352                                0.0, 1.0]; 
     353                        //options: L_theta = log value of theta, 
     354                        opt = "L_theta"; 
     355                }; 
     356                \endcode 
     357 
     358                Result is ARX data source offering with full history as Drv. 
     359                */ 
     360                void from_setting ( const Setting &set ); 
     361 
     362                // TODO dodelat void to_setting( Setting &set ) const; 
    267363}; 
    268364 
    … …  
    271367 
    272368class stateDS : public DS { 
    273 private: 
    274         //!conditional pdf of the state evolution \f$ f(x_t|x_{t-1}) \f$ 
    275         shared_ptr<mpdf> IM; 
    276  
    277         //!conditional pdf of the observations \f$ f(d_t|x_t) \f$ 
    278         shared_ptr<mpdf> OM; 
    279  
    280 protected: 
    281         //! result storage 
    282         vec dt; 
    283         //! state storage 
    284         vec xt; 
    285         //! input storage 
    286         vec ut; 
    287         //! Logger 
    288         int L_xt; 
    289  
    290 public: 
    291         void getdata ( vec &dt0 ) { 
    292                 dt0 = dt; 
    293         } 
    294  
    295         void getdata ( vec &dt0, const ivec &indices ) { 
    296                 dt0 = dt ( indices ); 
    297         } 
    298  
    299         stateDS ( const shared_ptr<mpdf> &IM0, const shared_ptr<mpdf> &OM0, int usize ) : IM ( IM0 ), OM ( OM0 ), 
    300                 dt ( OM0->dimension() ), xt ( IM0->dimension() ), 
    301                 ut ( usize ), L_xt(0) { } 
    302  
    303         stateDS() : L_xt(0) { } 
    304  
    305         virtual void step() { 
    306                 xt = IM->samplecond ( concat ( xt, ut ) ); 
    307                 dt = OM->samplecond ( concat ( xt, ut ) ); 
    308         } 
    309  
    310         virtual void log_add ( logger &L ) { 
    311                 DS::log_add ( L ); 
    312                 L_xt = L.add ( IM->_rv(), "true" ); 
    313         } 
    314         virtual void logit ( logger &L ) { 
    315                 DS::logit ( L ); 
    316                 L.logit ( L_xt, xt ); 
    317         } 
    318  
    319         /*! UI for stateDS 
    320  
    321         The DS is constructed from a structure with fields: 
    322         \code 
    323         system = { 
    324                 type = "stateDS"; 
    325                 //Internal model 
    326                 IM = { type = "mpdf"; //<-- valid offspring! e.g. "mlnorm" 
    327                         rv = { //description of x_t 
    328                                 names=["name1",...]; 
    329                                 sizes=[2,1]; // optional default=[1,1...]; 
    330                                 times=[0,0]; // optional default=[0,0...]; 
    331                                 } 
    332                         rvu= { //description of  u_t 
    333                                 //optional default=empty 
    334                                 } 
    335  
    336                         // remaining fields depending on the chosen type 
    337                         }; 
    338                 //Observation model 
    339                 OM = { type = "mpdf-offspring"; 
    340                         rv = {}; //description of d_t 
    341                         rvu = {type="internal", path="system.IM.rvu"}; //description of u_t 
    342  
    343                         //remaining fields 
    344                 } 
    345         }; 
    346         \endcode 
    347         */ 
    348         void from_setting ( const Setting &set ); 
    349  
    350         // TODO dodelat void to_setting( Setting &set ) const; 
     369        private: 
     370                //!conditional pdf of the state evolution \f$ f(x_t|x_{t-1}) \f$ 
     371                shared_ptr<mpdf> IM; 
     372 
     373                //!conditional pdf of the observations \f$ f(d_t|x_t) \f$ 
     374                shared_ptr<mpdf> OM; 
     375 
     376        protected: 
     377                //! result storage 
     378                vec dt; 
     379                //! state storage 
     380                vec xt; 
     381                //! input storage 
     382                vec ut; 
     383                //! Logger 
     384                int L_xt; 
     385 
     386        public: 
     387                void getdata ( vec &dt0 ) { 
     388                        dt0 = dt; 
     389                } 
     390 
     391                void getdata ( vec &dt0, const ivec &indices ) { 
     392                        dt0 = dt ( indices ); 
     393                } 
     394 
     395                stateDS ( const shared_ptr<mpdf> &IM0, const shared_ptr<mpdf> &OM0, int usize ) : IM ( IM0 ), OM ( OM0 ), 
     396                                dt ( OM0->dimension() ), xt ( IM0->dimension() ), 
     397                                ut ( usize ), L_xt ( 0 ) { } 
     398 
     399                stateDS() : L_xt ( 0 ) { } 
     400 
     401                virtual void step() { 
     402                        xt = IM->samplecond ( concat ( xt, ut ) ); 
     403                        dt = OM->samplecond ( concat ( xt, ut ) ); 
     404                } 
     405 
     406                virtual void log_add ( logger &L ) { 
     407                        DS::log_add ( L ); 
     408                        L_xt = L.add ( IM->_rv(), "true" ); 
     409                } 
     410                virtual void logit ( logger &L ) { 
     411                        DS::logit ( L ); 
     412                        L.logit ( L_xt, xt ); 
     413                } 
     414 
     415                /*! UI for stateDS 
     416 
     417                The DS is constructed from a structure with fields: 
     418                \code 
     419                system = { 
     420                        type = "stateDS"; 
     421                        //Internal model 
     422                        IM = { type = "mpdf"; //<-- valid offspring! e.g. "mlnorm" 
     423                                rv = { //description of x_t 
     424                                        names=["name1",...]; 
     425                                        sizes=[2,1]; // optional default=[1,1...]; 
     426                                        times=[0,0]; // optional default=[0,0...]; 
     427                                        } 
     428                                rvu= { //description of  u_t 
     429                                        //optional default=empty 
     430                                        } 
     431 
     432                                // remaining fields depending on the chosen type 
     433                                }; 
     434                        //Observation model 
     435                        OM = { type = "mpdf-offspring"; 
     436                                rv = {}; //description of d_t 
     437                                rvu = {type="internal", path="system.IM.rvu"}; //description of u_t 
     438 
     439                                //remaining fields 
     440                        } 
     441                }; 
     442                \endcode 
     443                */ 
     444                void from_setting ( const Setting &set ); 
     445 
     446                // TODO dodelat void to_setting( Setting &set ) const; 
    351447 
    352448}; 

  • library/bdm/design/ctrlbase.h

    r586 r598  
    9696                        pre_qr.set_submatrix(0,0,s*pr); 
    9797                        pre_qr.set_submatrix(dimx+dimu+dimy, dimu+dimx, -Qy*y_req); 
    98                         if (!qr(pre_qr,post_qr)) bdm_warning("QR in LQG unstable"); 
     98                        if (!qr(pre_qr,post_qr)){ bdm_warning("QR in LQG unstable");} 
    9999                        triu(post_qr); 
    100100        // hn(m+1:2*m+n+r,m+1:2*m+n+r); 

  • library/tests/CMakeLists.txt

    r586 r598  
    3939 
    4040# using UnitTest++ 
    41 add_executable(testsuite datalink_test.cpp egiw_test.cpp emix_test.cpp epdf_test.cpp logger_test.cpp LQG_test.cpp merger_test.cpp  
     41add_executable(testsuite datalink_test.cpp  datasource_test.cpp egiw_test.cpp emix_test.cpp epdf_test.cpp logger_test.cpp LQG_test.cpp merger_test.cpp  
    4242        mpdf_test.cpp randun_test.cpp rectangular_support_test.cpp rv_test.cpp shared_ptr_test.cpp square_mat_test.cpp testsuite.cpp user_info_test.cpp) 
    4343target_link_libraries(testsuite testutil unittest) 

  • library/tests/datalink_test.cpp

    r545 r598  
    8787        } 
    8888} 
     89 
     90TEST ( test_datalink_buffered ) { 
     91        RV a = RV ( "{dl_a }" ); 
     92        RV b = RV ( "{dlb }" ); 
     93        RV ab=concat(a,b); 
     94        RV aa = a; 
     95        a.t(-1); 
     96        aa.add(a);  
     97 
     98         
     99        datalink_buffered dl; 
     100        dl.set_connection ( concat(aa,b), ab); 
     101 
     102        vec val_up ( "1 2" ); 
     103        dl.step(val_up); 
     104        val_up="3 4"; 
     105         
     106         
     107        vec p = dl.pushdown ( val_up ); 
     108        vec exp_p = " 3, 1, 4 "; 
     109        CHECK_EQUAL ( exp_p, p ); 
     110 
     111} 

  • library/tests/epdf_harness.cpp

    r536 r598  
    8181                // test of pdflog at zero 
    8282                vec zero ( 0 ); 
    83                 vec zeron ( hepdf->dimension() ); 
    84                 for ( int i = 0; i < zeron.size(); ++i ) { 
    85                         zeron ( i ) = 0; 
    86                 } 
     83                vec zeron=zeros ( hepdf->dimension() ); 
    8784 
    8885                double lpz = hepdf->evallog ( zeron ); 
    89                 CHECK_CLOSE_EX ( lpz, mEp.evallogcond ( zeron, zero ), tolerance ); 
     86                double lpzc=mEp.evallogcond ( zeron, zero ); 
     87                CHECK_CLOSE_EX ( lpz, lpzc, tolerance ); 
    9088 
    9189                vec lpzv(1);