Changeset 271 for bdm

Show
Ignore:
Timestamp:
02/16/09 10:03:13 (15 years ago)
Author:
smidl
Message:

Next major revision

Location:
bdm
Files:
17 modified

Legend:

Unmodified
Added
Removed
  • bdm/estim/arx.h

    r270 r271  
    3131\f] 
    3232 
    33 Extension for time-variant parameters \f$\theta_t,r_t\f$ may be achived using exponential forgetting (Kulhavy and Zarrop, 1993). In such a case, the forgetting factor \c frg \f$\in <0,1>\f$ should be given in the constructor. Time-invariant parameters are estimated for \c frg = 1. 
     33See \ref tut_arx for mathematical treatment. 
     34 
     35The easiest way how to use the class is: 
     36\include arx_simple.cpp 
     37 
    3438*/ 
    3539class ARX: public BMEF { 
     
    8993        //!@{ 
    9094        const egiw* _e() const {return &est ;}; 
    91         const egiw& _epdf() const {return est;} 
     95        const egiw& posterior() const {return est;} 
    9296        //!@} 
    9397 
  • bdm/estim/libKF.cpp

    r270 r271  
    66using std::endl; 
    77 
    8 KalmanFull::KalmanFull ( mat A0, mat B0, mat C0, mat D0, mat R0, mat Q0, mat P0, vec mu0 ) { 
     8KalmanFull::KalmanFull ( mat A0, mat B0, mat C0, mat D0, mat Q0, mat R0, mat P0, vec mu0 ) { 
    99        dimx = A0.rows(); 
    1010        dimu = B0.cols(); 
     
    1515        it_assert_debug ( C0.cols() ==dimx, "KalmanFull: C is not square" ); 
    1616        it_assert_debug ( ( D0.rows() ==dimy ) || ( D0.cols() ==dimu ), "KalmanFull: D is not compatible" ); 
     17        it_assert_debug ( ( Q0.cols() ==dimx ) || ( Q0.rows() ==dimx ), "KalmanFull: Q is not compatible" ); 
    1718        it_assert_debug ( ( R0.cols() ==dimy ) || ( R0.rows() ==dimy ), "KalmanFull: R is not compatible" ); 
    18         it_assert_debug ( ( Q0.cols() ==dimx ) || ( Q0.rows() ==dimx ), "KalmanFull: Q is not compatible" ); 
    1919 
    2020        A = A0; 
     
    116116 
    117117 
    118 void KalmanCh::set_parameters ( const mat &A0,const mat &B0,const mat &C0,const mat &D0,const chmat &R0,const chmat &Q0 ) { 
    119  
    120         ( ( Kalman<chmat>* ) this )->set_parameters ( A0,B0,C0,D0,R0,Q0 ); 
     118void KalmanCh::set_parameters ( const mat &A0,const mat &B0,const mat &C0,const mat &D0,const chmat &Q0,const chmat &R0 ) { 
     119 
     120        ( ( Kalman<chmat>* ) this )->set_parameters ( A0,B0,C0,D0,Q0,R0 ); 
    121121        // Cholesky special! 
    122         preA.clear(); 
     122        preA=zeros(dimy+dimx+dimx,dimy+dimx); 
     123//      preA.clear(); 
    123124        preA.set_submatrix ( 0,0,R._Ch() ); 
    124125        preA.set_submatrix ( dimy+dimx,dimy,Q._Ch() ); 
  • bdm/estim/libKF.h

    r270 r271  
    109109        Kalman ( const Kalman<sq_T> &K0 ); 
    110110        //! Set parameters with check of relevance 
    111         void set_parameters ( const mat &A0,const mat &B0,const mat &C0,const mat &D0,const sq_T &R0,const sq_T &Q0 ); 
     111        void set_parameters ( const mat &A0,const mat &B0,const mat &C0,const mat &D0,const sq_T &Q0,const sq_T &R0 ); 
    112112        //! Set estimate values, used e.g. in initialization. 
    113113        void set_est ( const vec &mu0, const sq_T &P0 ) { 
     
    121121        void bayes ( const vec &dt ); 
    122122        //!access function 
    123         const epdf& _epdf() const {return est;} 
     123        const epdf& posterior() const {return est;} 
    124124        const enorm<sq_T>* _e() const {return &est;} 
    125125        //!access function 
     
    130130 
    131131/*! \brief Kalman filter in square root form 
     132 
     133Trivial example: 
     134\include kalman_simple.cpp 
     135  
    132136*/ 
    133137class KalmanCh : public Kalman<chmat>{ 
     
    142146        KalmanCh ():Kalman<chmat>(),preA(),postA(){}; 
    143147        //! Set parameters with check of relevance 
    144         void set_parameters ( const mat &A0,const mat &B0,const mat &C0,const mat &D0,const chmat &R0,const chmat &Q0 ); 
    145         void set_est ( const vec &mu0, const chmat &P0 ) { 
     148        void set_parameters ( const mat &A0,const mat &B0,const mat &C0,const mat &D0,const chmat &Q0,const chmat &R0 ); 
     149        void set_statistics ( const vec &mu0, const chmat &P0 ) { 
    146150                est.set_parameters ( mu0,P0 ); 
    147151        }; 
     
    187191        void set_est (vec mu0, mat P0){mu=mu0;P=P0;}; 
    188192        //!dummy! 
    189         const epdf& _epdf()const{return E;}; 
     193        const epdf& posterior()const{return E;}; 
    190194        const enorm<fsqmat>* _e()const{return &E;}; 
    191195        const mat _R(){return P;} 
     
    281285 
    282286template<class sq_T> 
    283 void Kalman<sq_T>::set_parameters ( const mat &A0,const  mat &B0, const mat &C0, const mat &D0, const sq_T &R0, const sq_T &Q0 ) { 
     287void Kalman<sq_T>::set_parameters ( const mat &A0,const  mat &B0, const mat &C0, const mat &D0, const sq_T &Q0, const sq_T &R0 ) { 
     288        dimx = A0.rows(); 
     289        dimy = C0.rows(); 
     290        dimu = B0.cols(); 
     291         
    284292        it_assert_debug ( A0.cols() ==dimx, "Kalman: A is not square" ); 
    285293        it_assert_debug ( B0.rows() ==dimx, "Kalman: B is not compatible" ); 
  • bdm/estim/libPF.h

    r270 r271  
    122122                for ( int i=0;i<n;i++ ) { 
    123123                        Bms[i] = new BM_T ( BMcond0 ); //copy constructor 
    124                         const epdf& pom=Bms[i]->_epdf(); 
     124                        const epdf& pom=Bms[i]->posterior(); 
    125125                        jest.set_elements ( i,1.0/n,&pom ); 
    126126                } 
     
    131131 
    132132        void bayes ( const vec &dt ); 
    133         const epdf& _epdf() const {return jest;} 
     133        const epdf& posterior() const {return jest;} 
    134134        const epdf* _e() const {return &jest;} //Fixme: is it useful? 
    135135        //! Set postrior of \c rvc to samples from epdf0. Statistics of Bms are not re-computed! Use only for initialization! 
     
    195195                                delete Bms[i]; 
    196196                                Bms[i] = new BM_T ( *Bms[ind ( i ) ] ); //copy constructor 
    197                                 const epdf& pom=Bms[i]->_epdf(); 
     197                                const epdf& pom=Bms[i]->posterior(); 
    198198                                jest.set_elements ( i,1.0/n,&pom ); 
    199199                        } 
  • bdm/estim/merger.h

    r270 r271  
    8181        void merge ( const epdf* g0 ); 
    8282//!Create a mixture density, make sure to call init() before the first call 
    83         void merge () {merge ( & ( Mix._epdf() ) );}; 
     83        void merge () {merge ( & ( Mix.posterior() ) );}; 
    8484 
    8585//! Merge log-likelihood values 
     
    8787//! sample from merged density 
    8888//! weight w is a 
    89         vec sample ( ) const { return Mix._epdf().sample();} 
     89        vec sample ( ) const { return Mix.posterior().sample();} 
    9090        double evallog ( const vec &dt ) const { 
    9191                vec dtf=ones ( dt.length() +1 ); 
  • bdm/estim/mixef.cpp

    r270 r271  
    124124double MixEF::logpred ( const vec &dt ) const { 
    125125 
    126         vec w=weights._epdf().mean(); 
     126        vec w=weights.posterior().mean(); 
    127127        double exLL=0.0; 
    128128        for ( int i=0;i<n;i++ ) { 
     
    137137        emix* tmp; 
    138138        tmp = new emix (  ); 
    139         tmp->set_parameters ( weights._epdf().mean(), pC, false ); 
     139        tmp->set_parameters ( weights.posterior().mean(), pC, false ); 
    140140        tmp->ownComs(); 
    141141        return tmp; 
  • bdm/estim/mixef.h

    r270 r271  
    5959                for ( int i=0;i<Coms.length();i++ ) { 
    6060//                      it_assert_debug(!x,"MixEF::MixEF : Incompatible components"); 
    61                         epdfs ( i ) =& ( Coms ( i )->_epdf() ); 
     61                        epdfs ( i ) =& ( Coms ( i )->posterior() ); 
    6262                } 
    6363                // last in the product is the weight 
    64                 epdfs ( n ) =& ( weights._epdf() ); 
     64                epdfs ( n ) =& ( weights.posterior() ); 
    6565                est = new eprod ( epdfs ); 
    6666        } 
     
    104104        void bayesB ( const mat &dt, const vec &wData ); 
    105105        double logpred ( const vec &dt ) const; 
    106         const epdf& _epdf() const {return *est;} 
     106        const epdf& posterior() const {return *est;} 
    107107        const eprod* _e() const {return est;} 
    108108        emix* epredictor() const; 
  • bdm/stat/emix.cpp

    r270 r271  
    6060//                      // If rvaddok==false, mpdfs overlap => assert error. 
    6161//                      it_assert_debug(rvaddok||overlap,"mprod::mprod() input mpdfs overlap in rv!"); 
    62 //                      epdfs ( i ) = & ( mpdfs ( i )->_epdf() ); // add pointer to epdf 
     62//                      epdfs ( i ) = & ( mpdfs ( i )->posterior() ); // add pointer to epdf 
    6363//              }; 
    6464//              // Create rvc 
  • bdm/stat/emix.h

    r270 r271  
    163163class mprod: public compositepdf, public mpdf { 
    164164protected: 
    165         //! pointers to epdfs - shortcut to mpdfs()._epdf() 
     165        //! pointers to epdfs - shortcut to mpdfs().posterior() 
    166166        Array<epdf*> epdfs; 
    167167        //! Data link for each mpdfs 
  • bdm/stat/libBM.cpp

    r270 r271  
    55//! Space of basic BDM structures 
    66namespace bdm { 
    7          
     7 
    88const int RV_BUFFER_STEP=1; 
    99RVmap RV_MAP; 
    10 Array<string> RV_NAMES(RV_BUFFER_STEP); 
    11 ivec RV_SIZES(RV_BUFFER_STEP); 
     10Array<string> RV_NAMES ( RV_BUFFER_STEP ); 
     11ivec RV_SIZES ( RV_BUFFER_STEP ); 
    1212 
    1313RV RV0=RV(); 
     
    1818        RVmap::const_iterator iter = RV_MAP.find ( name ); 
    1919        if ( iter == RV_MAP.end() ) { 
    20                 id=RV_MAP.size()+1; 
     20                id=RV_MAP.size() +1; 
    2121                RV_MAP.insert ( make_pair ( name,id ) ); //add new rv 
    22                 if (id>=RV_NAMES.length()){ 
    23                         RV_NAMES.set_length(id+RV_BUFFER_STEP,true); 
    24                         RV_SIZES.set_length(id+RV_BUFFER_STEP,true); 
    25                 } 
    26                 RV_NAMES(id)=name; 
    27                 RV_SIZES(id)=size; 
     22                if ( id>=RV_NAMES.length() ) { 
     23                        RV_NAMES.set_length ( id+RV_BUFFER_STEP,true ); 
     24                        RV_SIZES.set_length ( id+RV_BUFFER_STEP,true ); 
     25                } 
     26                RV_NAMES ( id ) =name; 
     27                RV_SIZES ( id ) =size; 
    2828        } 
    2929        else { 
    3030                id = iter->second; 
    31                 it_warning ("RV of name " + name + "already exists" ); 
     31                it_assert(RV_SIZES(id)==size,"RV "+name+" of different size already exists"); 
    3232        } 
    3333        return id; 
    3434}; 
    3535 
    36 int RV::countsize() const{ int tmp=0; for(int i=0;i<len;i++){tmp+=RV_SIZES(ids(i));} return tmp;} 
     36int RV::countsize() const{ int tmp=0; for ( int i=0;i<len;i++ ) {tmp+=RV_SIZES ( ids ( i ) );} return tmp;} 
     37 
     38ivec RV::cumsizes() const { 
     39        ivec szs ( len ); 
     40        int tmp=0; 
     41        for ( int i=0;i<len;i++ ) { 
     42                tmp+=RV_SIZES ( ids ( i ) ); 
     43                szs ( i ) = tmp; 
     44        } 
     45        return szs; 
     46} 
    3747 
    3848void RV::init ( Array<std::string> in_names, ivec in_sizes,ivec in_times ) { 
     
    4050        it_assert_debug ( in_names.length() ==in_times.length(), "check \"times\" " ); 
    4151        it_assert_debug ( in_names.length() ==in_sizes.length(), "check \"sizes\" " ); 
    42          
     52 
    4353        times.set_length ( len ); 
    4454        ids.set_length ( len ); 
     
    5464RV::RV ( string name, int sz, int tm ) { 
    5565        Array<string> A ( 1 ); A ( 0 ) =name; 
    56         init (A,vec_1 ( sz ),vec_1 ( tm ) ); 
     66        init ( A,vec_1 ( sz ),vec_1 ( tm ) ); 
    5767} 
    5868 
     
    8292RV RV::subselect ( const ivec &ind ) const { 
    8393        RV ret; 
    84         ret.ids = ids(ind); 
    85         ret.times= times(ind); 
     94        ret.ids = ids ( ind ); 
     95        ret.times= times ( ind ); 
    8696        ret.len = ind.length(); 
    8797        ret.dsize=ret.countsize(); 
     
    118128        int pos = 0; 
    119129        for ( i = 0;i < len;i++ ) { 
    120                 idlist.set_subvector ( pos, pos + size ( ids(i) ) - 1, ids ( i ) ); 
    121                 tmlist.set_subvector ( pos, pos + size ( ids(i) ) - 1, times ( i ) ); 
    122                 pos += size ( ids(i) ); 
     130                idlist.set_subvector ( pos, pos + size ( ids ( i ) ) - 1, ids ( i ) ); 
     131                tmlist.set_subvector ( pos, pos + size ( ids ( i ) ) - 1, times ( i ) ); 
     132                pos += size ( ids ( i ) ); 
    123133        } 
    124134        return str ( idlist, tmlist ); 
  • bdm/stat/libBM.h

    r270 r271  
    1111*/ 
    1212 
    13 /*! 
    14 \defgroup core Core BDM classes 
    15 @{ 
    16 \defgroup const Constructors 
    17 \defgroup math Mathematical operations 
    18 */ 
    1913#ifndef BM_H 
    2014#define BM_H 
     
    106100        int init ( const  string &name, int size ); 
    107101public: 
    108         //! \name Constructors  
    109         //!@{ 
    110          
    111         //! Full constructor  
     102        //! \name Constructors 
     103        //!@{ 
     104 
     105        //! Full constructor 
    112106        RV ( Array<std::string> in_names, ivec in_sizes, ivec in_times ) {init ( in_names,in_sizes,in_times );}; 
    113         //! Constructor with times=0  
     107        //! Constructor with times=0 
    114108        RV ( Array<std::string> in_names, ivec in_sizes ) {init ( in_names,in_sizes,zeros_i ( in_names.length() ) );}; 
    115         //! Constructor with sizes=1, times=0  
     109        //! Constructor with sizes=1, times=0 
    116110        RV ( Array<std::string> in_names ) {init ( in_names,ones_i ( in_names.length() ),zeros_i ( in_names.length() ) );} 
    117         //! Constructor of empty RV  
     111        //! Constructor of empty RV 
    118112        RV () :dsize ( 0 ),len ( 0 ),ids ( 0 ),times ( 0 ) {}; 
    119113        //! Constructor of a single RV with given id 
    120114        RV ( string name, int sz, int tm=0 ); 
    121115        //!@} 
    122          
     116 
    123117        //! \name Access functions 
    124118        //!@{ 
    125          
     119 
    126120        //! Printing output e.g. for debugging. 
    127121        friend std::ostream &operator<< ( std::ostream &os, const RV &rv ); 
     
    129123        //! Recount size of the corresponding data vector 
    130124        int countsize() const; 
     125        ivec cumsizes() const; 
    131126        int length() const {return len;} ; 
    132127        int id ( int at ) const{return ids ( at );}; 
    133         int size ( int at ) const {return RV_SIZES ( at );}; 
     128        int size ( int at ) const {return RV_SIZES ( ids(at) );}; 
    134129        int time ( int at ) const{return times ( at );}; 
    135         std::string name ( int at ) const {return RV_NAMES ( at );}; 
     130        std::string name ( int at ) const {return RV_NAMES ( ids(at) );}; 
    136131        void set_time ( int at, int time0 ) {times ( at ) =time0;}; 
    137132        //!@} 
    138          
     133 
    139134        //TODO why not inline and later?? 
    140135 
    141136        //! \name Algebra on Random Variables 
    142137        //!@{ 
    143          
     138 
    144139        //! Find indices of self in another rv, \return ivec of the same size as self. 
    145140        ivec findself ( const RV &rv2 ) const; 
     
    154149        //! Select only variables at indeces ind 
    155150        RV operator() ( const ivec &ind ) const {return subselect ( ind );}; 
     151        //! Select from data vector starting at di1 to di2 
     152        RV operator() ( int di1, int di2 ) const { 
     153                ivec sz=cumsizes(); 
     154                int i1=0; 
     155                while ( sz ( i1 ) <di1 ) i1++; 
     156                int i2=i1; 
     157                while ( sz ( i2 ) <di2 ) i2++; 
     158                return subselect ( linspace ( i1,i2 ) ); 
     159        }; 
    156160        //! Shift \c time shifted by delta. 
    157161        void t ( int delta ); 
    158162        //!@} 
    159          
    160         //!\name Relation to vectors  
    161         //!@{ 
    162          
     163 
     164        //!\name Relation to vectors 
     165        //!@{ 
     166 
    163167        //! generate \c str from rv, by expanding sizes 
    164168        str tostr() const; 
     
    172176        int mint () const {return min ( times );}; 
    173177        //!@} 
    174          
     178 
    175179}; 
    176180 
     
    216220public: 
    217221        /*! \name Constructors 
    218          Construction of each epdf should support two types of constructors:  
    219         \li empty constructor,  
     222         Construction of each epdf should support two types of constructors: 
     223        \li empty constructor, 
    220224        \li copy constructor, 
    221          
     225 
    222226        The following constructors should be supported for convenience: 
    223         \li constructor followed by calling \c set_parameters()  
     227        \li constructor followed by calling \c set_parameters() 
    224228        \li constructor accepting random variables calling \c set_rv() 
    225          
     229 
    226230         All internal data structures are constructed as empty. Their values (including sizes) will be set by method \c set_parameters(). This way references can be initialized in constructors. 
    227231        @{*/ 
    228         epdf() :dim(0),rv ( ) {}; 
    229         epdf(const epdf &e) :dim(e.dim),rv (e.rv) {}; 
    230         epdf(const RV &rv0) {set_rv(rv0);}; 
    231         void set_parameters(int dim0){dim=dim0;} 
    232         //!@} 
    233          
     232        epdf() :dim ( 0 ),rv ( ) {}; 
     233        epdf ( const epdf &e ) :dim ( e.dim ),rv ( e.rv ) {}; 
     234        epdf ( const RV &rv0 ) {set_rv ( rv0 );}; 
     235        void set_parameters ( int dim0 ) {dim=dim0;} 
     236        //!@} 
     237 
    234238        //! \name Matematical Operations 
    235239        //!@{ 
    236          
     240 
    237241        //! Returns a sample, \f$ x \f$ from density \f$ f_x()\f$ 
    238         virtual vec sample () const {it_error("not implemneted");return vec(0);}; 
     242        virtual vec sample () const {it_error ( "not implemneted" );return vec ( 0 );}; 
    239243        //! Returns N samples, \f$ [x_1 , x_2 , \ldots \ \f$  from density \f$ f_x(rv)\f$ 
    240244        virtual mat sample_m ( int N ) const; 
    241245        //! Compute log-probability of argument \c val 
    242         virtual double evallog ( const vec &val ) const {it_error("not implemneted");return 0.0;}; 
     246        virtual double evallog ( const vec &val ) const {it_error ( "not implemneted" );return 0.0;}; 
    243247        //! Compute log-probability of multiple values argument \c val 
    244248        virtual vec evallog_m ( const mat &Val ) const { 
     
    252256        virtual epdf* marginal ( const RV &rv ) const {it_warning ( "Not implemented" ); return NULL;} 
    253257        //! return expected value 
    254         virtual vec mean() const {it_error("not implemneted");return vec(0);}; 
     258        virtual vec mean() const {it_error ( "not implemneted" );return vec ( 0 );}; 
    255259        //! return expected variance (not covariance!) 
    256         virtual vec variance() const {it_error("not implemneted");return vec(0);}; 
    257         //!@} 
    258          
     260        virtual vec variance() const {it_error ( "not implemneted" );return vec ( 0 );}; 
     261        //!@} 
     262 
    259263        //! \name Connection to other classes 
    260         //! Description of the random quantity via attribute \c rv is optional.  
    261         //! For operations such as sampling \c rv does not need to be set. However, for \c marginalization  
    262         //! and \c conditioning \c rv has to be set. NB:  
     264        //! Description of the random quantity via attribute \c rv is optional. 
     265        //! For operations such as sampling \c rv does not need to be set. However, for \c marginalization 
     266        //! and \c conditioning \c rv has to be set. NB: 
    263267        //! @{ 
    264          
     268 
    265269        //!Name its rv 
    266270        void set_rv ( const RV &rv0 ) {rv = rv0; }//it_assert_debug(isnamed(),""); }; 
    267         //! True if rv is assigned  
    268         bool isnamed() const {bool b=( dim==rv._dsize() );return b;} 
     271        //! True if rv is assigned 
     272        bool isnamed() const {bool b= ( dim==rv._dsize() );return b;} 
    269273        //! Return name (fails when isnamed is false) 
    270274        const RV& _rv() const {it_assert_debug ( isnamed(),"" ); return rv;} 
    271275        //!@} 
    272          
     276 
    273277        //! \name Access to attributes 
    274278        //! @{ 
    275          
     279 
    276280        //! Size of the random variable 
    277281        int dimension() const {return dim;} 
    278282        //!@} 
    279          
     283 
    280284}; 
    281285 
     
    295299        //! \name Constructors 
    296300        //! @{ 
    297          
    298         mpdf ( ) :dimc(0),rvc ( ) {}; 
     301 
     302        mpdf ( ) :dimc ( 0 ),rvc ( ) {}; 
    299303        //! copy constructor does not set pointer \c ep - has to be done in offsprings! 
    300         mpdf (const mpdf &m ) :dimc(m.dimc),rvc (m.rvc ) {}; 
     304        mpdf ( const mpdf &m ) :dimc ( m.dimc ),rvc ( m.rvc ) {}; 
    301305        //!@} 
    302306 
    303307        //! \name Matematical operations 
    304308        //!@{ 
    305          
     309 
    306310        //! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv 
    307311        virtual vec samplecond ( const vec &cond ) { 
     
    330334        //! \name Access to attributes 
    331335        //! @{ 
    332          
     336 
    333337        RV _rv() {return ep->_rv();} 
    334338        RV _rvc() {it_assert_debug ( isnamed(),"" ); return rvc;} 
     
    338342        epdf* _e() {return ep;} 
    339343        //!@} 
    340          
     344 
    341345        //! \name Connection to other objects 
    342346        //!@{ 
    343347        void set_rvc ( const RV &rvc0 ) {rvc=rvc0;} 
    344         void set_rv ( const RV &rv0 ) {ep->set_rv(rv0);} 
    345         bool isnamed() {return (ep->isnamed())&&(dimc==rvc._dsize());} 
     348        void set_rv ( const RV &rv0 ) {ep->set_rv ( rv0 );} 
     349        bool isnamed() {return ( ep->isnamed() ) && ( dimc==rvc._dsize() );} 
    346350        //!@} 
    347351}; 
     
    382386public: 
    383387        //! Constructor 
    384         datalink ( const RV &rv, const RV &rv_up ) : 
    385                         downsize ( rv._dsize() ), upsize ( rv_up._dsize() ), v2v_up ( rv.dataind ( rv_up ) )  { 
     388        datalink () {}; 
     389        datalink ( const RV &rv, const RV &rv_up ) {set_connection ( rv,rv_up );}; 
     390        //! set connection, rv must be fully present in rv_up 
     391        void set_connection ( const RV &rv, const RV &rv_up ) { 
     392                downsize = rv._dsize(); 
     393                upsize = rv_up._dsize(); 
     394                v2v_up= ( rv.dataind ( rv_up ) ); 
     395 
    386396                it_assert_debug ( v2v_up.length() ==downsize,"rv is not fully in rv_up" ); 
    387397        } 
     
    471481        logger ( ) : entries ( 0 ),names ( 0 ) {} 
    472482 
    473         //! returns an identifier which will be later needed for calling the log() function 
     483        //! returns an identifier which will be later needed for calling the \c logit() function 
     484        //! For empty RV it returns -1, this entry will be ignored by \c logit(). 
    474485        virtual int add ( const RV &rv, string name="" ) { 
    475                 int id=entries.length(); 
    476                 names=concat ( names, name ); // diff 
    477                 entries.set_length ( id+1,true ); 
    478                 entries ( id ) = rv; 
     486                int id; 
     487                if ( rv._dsize() >0 ) { 
     488                        id=entries.length(); 
     489                        names=concat ( names, name ); // diff 
     490                        entries.set_length ( id+1,true ); 
     491                        entries ( id ) = rv; 
     492                } 
     493                else { id =-1;} 
    479494                return id; // identifier of the last entry 
    480495        } 
     
    539554public: 
    540555        //! default constructors 
    541         DS() :Drv (  ),Urv ( ) {}; 
     556        DS() :Drv ( ),Urv ( ) {}; 
    542557        //! Returns full vector of observed data=[output, input] 
    543558        virtual void getdata ( vec &dt ) {it_error ( "abstract class" );}; 
     
    590605        //! \name Constructors 
    591606        //! @{ 
    592          
    593         BM () :ll (0),evalll ( false) {}; 
    594         BM ( const BM &B ) :  drv(B.drv), ll ( B.ll ), evalll ( B.evalll ) {} 
     607 
     608        BM () :ll ( 0 ),evalll ( false ) {}; 
     609        BM ( const BM &B ) :  drv ( B.drv ), ll ( B.ll ), evalll ( B.evalll ) {} 
    595610        //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually! 
    596611        //! Prototype: \code BM* _copy_(){return new BM(*this);} \endcode 
     
    600615        //! \name Mathematical operations 
    601616        //!@{ 
    602          
     617 
    603618        /*! \brief Incremental Bayes rule 
    604619        @param dt vector of input data 
     
    614629 
    615630        //!Constructs a predictive density \f$ f(d_{t+1} |d_{t}, \ldots d_{0}) \f$ 
    616         virtual epdf* epredictor (  ) const {it_error ( "Not implemented" );return NULL;}; 
     631        virtual epdf* epredictor ( ) const {it_error ( "Not implemented" );return NULL;}; 
    617632        //!Constructs a conditional density 1-step ahead predictor \f$ f(d_{t+1} |d_{t+h-1}, \ldots d_{t}) 
    618633        virtual mpdf* predictor ( ) const {it_error ( "Not implemented" );return NULL;}; 
    619634        //!@} 
    620          
     635 
    621636        //! \name Access to attributes 
    622637        //!@{ 
    623          
     638 
    624639        const RV& _drv() const {return drv;} 
    625640        void set_drv ( const RV &rv ) {drv=rv;} 
     641        void set_rv ( const RV &rv ) {const_cast<epdf&> ( posterior() ).set_rv ( rv );} 
    626642        double _ll() const {return ll;} 
    627643        void set_evalll ( bool evl0 ) {evalll=evl0;} 
    628         virtual const epdf& _epdf() const =0; 
     644        virtual const epdf& posterior() const =0; 
    629645        virtual const epdf* _e() const =0; 
    630646        //!@} 
     
    657673 
    658674}; //namespace 
    659 /*! @} */ 
    660675#endif // BM_H 
  • bdm/stat/libDS.cpp

    r270 r271  
    4343void ArxDS::step() { 
    4444        //shift history 
    45         H.shift_right ( 0, Drv._dsize() +Urv._dsize() ); 
     45        H.shift_right ( 0, dt_size ); 
    4646 
    47         H.set_subvector ( Drv._dsize(), U ); // write U after Drv 
     47        H.set_subvector ( dt_size-utsize, U ); // write U after Drv 
    4848 
    4949        //get regressor 
    50         rgr = rgrlnk->pushdown ( H ); 
     50        rgr = rgrlnk.pushdown ( H ); 
    5151        // Eval Y 
    5252        H.set_subvector ( 0, model.samplecond ( rgr ) ); 
  • bdm/stat/libDS.h

    r270 r271  
    2020 
    2121namespace bdm { 
    22         /*! 
    23         * \brief Memory storage of off-line data column-wise 
     22/*! 
     23* \brief Memory storage of off-line data column-wise 
    2424 
    25         The data are stored in an internal matrix \c Data . Each column of Data corresponds to one discrete time observation \f$t\f$. Access to this matrix is via indices \c rowid and \c delays. 
     25The data are stored in an internal matrix \c Data . Each column of Data corresponds to one discrete time observation \f$t\f$. Access to this matrix is via indices \c rowid and \c delays. 
    2626 
    27         The data can be loaded from a file. 
    28         */ 
    29         class MemDS : public DS { 
    30                 //! internal matrix of data 
    31                 mat Data; 
    32                 //! active column in the Data matrix 
    33                 int time; 
    34                 //!  vector of rows that are presented in Dt 
    35                 ivec rowid; 
    36                 //! vector of delays that are presented in Dt 
    37                 ivec delays; 
     27The data can be loaded from a file. 
     28*/ 
     29class MemDS : public DS { 
     30        //! internal matrix of data 
     31        mat Data; 
     32        //! active column in the Data matrix 
     33        int time; 
     34        //!  vector of rows that are presented in Dt 
     35        ivec rowid; 
     36        //! vector of delays that are presented in Dt 
     37        ivec delays; 
    3838 
    39         public: 
    40                 void getdata ( vec &dt ); 
    41                 void getdata ( vec &dt, const ivec &indeces ); 
    42                 void set_rvs ( RV &drv, RV &urv ); 
    43                 void write ( vec &ut ) {it_error ( "MemDS::write is not supported" );} 
    44                 void write ( vec &ut,ivec &indices ) {it_error ( "MemDS::write is not supported" );} 
    45                 void step(); 
    46                 //!Default constructor 
    47                 MemDS ( mat &Dat, ivec &rowid, ivec &delays ); 
     39public: 
     40        void getdata ( vec &dt ); 
     41        void getdata ( vec &dt, const ivec &indeces ); 
     42        void set_rvs ( RV &drv, RV &urv ); 
     43        void write ( vec &ut ) {it_error ( "MemDS::write is not supported" );} 
     44        void write ( vec &ut,ivec &indices ) {it_error ( "MemDS::write is not supported" );} 
     45        void step(); 
     46        //!Default constructor 
     47        MemDS ( mat &Dat, ivec &rowid, ivec &delays ); 
     48}; 
     49 
     50/*! 
     51\brief Generator of ARX data 
     52 
     53*/ 
     54class ArxDS : public DS { 
     55protected: 
     56        //! Rv of the regressor 
     57        RV Rrv; 
     58        //! History, ordered as \f$[y_t, u_t, y_{t-1 }, u_{t-1}, \ldots]\f$ 
     59        vec H; 
     60        //! (future) input 
     61        vec U; 
     62        //! temporary variable for regressor 
     63        vec rgr; 
     64        //! data link: H -> rgr 
     65        datalink rgrlnk; 
     66        //! model of Y - linear Gaussian 
     67        mlnorm<chmat> model; 
     68        //! options 
     69        bool opt_L_theta; 
     70        //! loggers 
     71        int L_theta; 
     72        int L_R; 
     73        int dt_size; 
     74public: 
     75        void getdata ( vec &dt ) { 
     76                //it_assert_debug ( dt.length() ==Drv.count(),"ArxDS" ); 
     77                dt=H; 
     78        }; 
     79        void getdata ( vec &dt, const ivec &indices ) { 
     80                it_assert_debug ( dt.length() ==indices.length(),"ArxDS" ); 
     81                dt=H ( indices ); 
     82        }; 
     83        void write ( vec &ut ) { 
     84                //it_assert_debug ( ut.length() ==Urv.count(),"ArxDS" ); 
     85                U=ut; 
     86        }; 
     87        void write ( vec &ut, const ivec &indices ) { 
     88                it_assert_debug ( ut.length() ==indices.length(),"ArxDS" ); 
     89                set_subvector ( U, indices,ut ); 
     90        }; 
     91        void step(); 
     92        //!Default constructor 
     93        ArxDS ( ) {}; 
     94        //! Set parameters of the internal model, H is maximum time delay 
     95        void set_parameters ( const mat &Th0, const vec mu0, const chmat &sqR0 ) 
     96        { model.set_parameters ( Th0, mu0, sqR0 );}; 
     97        //! Set 
     98        void set_drv(RV &yrv, RV &urv, RV &rrv){ 
     99                Rrv = rrv; 
     100                Urv = urv; 
     101                dt_size = yrv._dsize()+urv._dsize(); 
     102                 
     103                RV drv = concat(yrv,urv); 
     104                Drv = drv; 
     105                int td = rrv.mint(); 
     106                H.set_size(drv._dsize()*(-td+1)); 
     107                U.set_size(Urv._dsize()); 
     108                for (int i=-1;i>=td;i--){ 
     109                        drv.t(-1); 
     110                        Drv.add(drv); //shift u1 
     111                } 
     112                rgrlnk.set_connection(rrv,Drv); 
     113                 
     114                dtsize = Drv._dsize(); 
     115                utsize = Urv._dsize(); 
     116        } 
     117        //! set options from a string 
     118        void set_options ( const string &s ) { 
     119                opt_L_theta= ( s.find ( "L_theta" ) !=string::npos ); 
     120        }; 
     121        virtual void log_add ( logger &L ) { 
     122                //DS::log_add ( L ); too long!! 
     123                L_dt=L.add ( Drv(0,dt_size),"" ); 
     124                L_ut=L.add ( Urv,"" ); 
     125 
     126                mat &A =model._A(); 
     127                mat R =model._R(); 
     128                if ( opt_L_theta ) {L_theta=L.add ( RV ( "{th }", vec_1 ( A.rows() *A.cols() ) ),"t" );} 
     129                if ( opt_L_theta ) {L_R=L.add ( RV ( "{R }", vec_1 ( R.rows() *R.cols() ) ),"r" );} 
     130        } 
     131        virtual void logit ( logger &L ) { 
     132                //DS::logit ( L ); 
     133                L.logit( L_dt, H.left(dt_size)); 
     134                L.logit(L_ut, U); 
     135                 
     136                mat &A =model._A(); 
     137                mat R =model._R(); 
     138                if ( opt_L_theta ) {L.logit ( L_theta,vec ( A._data(), A.rows() *A.cols() ) );}; 
     139                if ( opt_L_theta ) {L.logit ( L_R, vec ( R._data(), R.rows() *R.rows() ) );}; 
     140        } 
     141 
     142}; 
     143 
     144class stateDS : public DS { 
     145protected: 
     146        //!conditional pdf of the state evolution \f$ f(x_t|x_{t-1}) \f$ 
     147        mpdf* IM; 
     148        //!conditional pdf of the observations \f$ f(d_t|x_t) \f$ 
     149        mpdf* OM; 
     150        //! result storage 
     151        vec dt; 
     152        //! state storage 
     153        vec xt; 
     154        //! input storage 
     155        vec ut; 
     156        //! Logger 
     157        int L_xt; 
     158public: 
     159        void getdata ( vec &dt0 ) {dt0=dt;} 
     160        void getdata ( vec &dt0, const ivec &indeces ) {dt0=dt ( indeces );} 
     161 
     162        stateDS ( mpdf* IM0, mpdf* OM0, int usize ) :DS ( ),IM ( IM0 ),OM ( OM0 ), 
     163                        dt ( OM0->dimension() ), xt ( IM0->dimension() ), ut ( usize ) {} 
     164        ~stateDS() {delete IM; delete OM;} 
     165        virtual void step() { 
     166                xt=IM->samplecond ( concat ( xt,ut ) ); 
     167                dt=OM->samplecond ( concat ( xt,ut ) ); 
    48168        }; 
    49169 
    50         /*! 
    51         \brief Generator of ARX data 
     170        virtual void log_add ( logger &L ) { 
     171                DS::log_add ( L ); 
     172                L_xt=L.add ( IM->_rv(),"true" ); 
     173        } 
     174        virtual void logit ( logger &L ) { 
     175                DS::logit ( L ); 
     176                L.logit ( L_xt,xt ); 
     177        } 
    52178 
    53         */ 
    54         class ArxDS : public DS { 
    55         protected: 
    56                 //! Rv of the regressor 
    57                 RV Rrv; 
    58                 //! Rv of the history (full regressor) 
    59                 RV Hrv; 
    60                 //! History, ordered as \f$[y_t, u_t, y_{t-1 }, u_{t-1}, \ldots]\f$ 
    61                 vec H; 
    62                 //! (future) input 
    63                 vec U; 
    64                 //! temporary variable for regressor 
    65                 vec rgr; 
    66                 //! data link: H -> rgr 
    67                 datalink *rgrlnk; 
    68                 //! model of Y - linear Gaussian 
    69                 mlnorm<chmat> model; 
    70                 //! options 
    71                 bool opt_L_theta; 
    72                 //! loggers 
    73                 int L_theta; 
    74                 int L_R; 
    75         public: 
    76                 void getdata ( vec &dt ) { 
    77                         //it_assert_debug ( dt.length() ==Drv.count(),"ArxDS" ); 
    78                         dt=H.left ( Urv._dsize() +Drv._dsize() ); 
    79                 }; 
    80                 void getdata ( vec &dt, const ivec &indices ) { 
    81                         it_assert_debug ( dt.length() ==indices.length(),"ArxDS" ); 
    82                         dt=H ( indices ); 
    83                 }; 
    84                 void write ( vec &ut ) { 
    85                         //it_assert_debug ( ut.length() ==Urv.count(),"ArxDS" ); 
    86                         U=ut; 
    87                 }; 
    88                 void write ( vec &ut, const ivec &indices ) { 
    89                         it_assert_debug ( ut.length() ==indices.length(),"ArxDS" ); 
    90                         set_subvector ( U, indices,ut ); 
    91                 }; 
    92                 void step(); 
    93                 //!Default constructor 
    94                 ArxDS ( ){}; 
    95                 //! Set parameters of the internal model 
    96                 void set_parameters ( const mat &Th0, const vec mu0, const chmat &sqR0 ) 
    97                 { model.set_parameters ( Th0, mu0, sqR0 ); }; 
    98                 //! set options from a string 
    99                 void set_options ( const string &s ) { 
    100                         opt_L_theta= ( s.find ( "L_theta" ) !=string::npos ); 
    101                 }; 
    102                 virtual void log_add ( logger &L ) { 
    103                         DS::log_add ( L ); 
    104                         mat &A =model._A(); 
    105                         mat R =model._R(); 
    106                         if ( opt_L_theta ) {L_theta=L.add ( RV ( "{theta }", vec_1 ( A.rows() *A.cols() ) ),"t" );} 
    107                         if ( opt_L_theta ) {L_R=L.add ( RV ( "{R }", vec_1 ( R.rows() *R.cols() ) ),"r" );} 
    108                 } 
    109                 virtual void logit ( logger &L ) { 
    110                         DS::logit ( L ); 
    111                         mat &A =model._A(); 
    112                         mat R =model._R(); 
    113                         if ( opt_L_theta ) {L.logit ( L_theta,vec ( A._data(), A.rows() *A.cols() ) );}; 
    114                         if ( opt_L_theta ) {L.logit ( L_R, vec ( R._data(), R.rows() *R.rows() ) );}; 
    115                 } 
    116  
    117         }; 
    118  
    119         class ARXDS : public ArxDS { 
    120         public: 
    121                 ARXDS ( ) : ArxDS ( ) {} 
    122  
    123                 void getdata ( vec &dt ) {dt=H;} 
    124                 void getdata ( vec &dt, const ivec &indeces ) {dt=H ( indeces );} 
    125                 virtual RV _drv() const {return Hrv;} 
    126  
    127         }; 
    128  
    129         class stateDS : public DS { 
    130         protected: 
    131                 //!conditional pdf of the state evolution \f$ f(x_t|x_{t-1}) \f$ 
    132                 mpdf* IM; 
    133                 //!conditional pdf of the observations \f$ f(d_t|x_t) \f$ 
    134                 mpdf* OM; 
    135                 //! result storage 
    136                 vec dt; 
    137                 //! state storage 
    138                 vec xt; 
    139                 //! input storage 
    140                 vec ut; 
    141                 //! Logger 
    142                 int L_xt; 
    143         public: 
    144                 void getdata ( vec &dt0 ) {dt0=dt;} 
    145                 void getdata ( vec &dt0, const ivec &indeces ) {dt0=dt ( indeces );} 
    146  
    147                 stateDS ( mpdf* IM0, mpdf* OM0, int usize ) :DS ( ),IM ( IM0 ),OM ( OM0 ), 
    148                                 dt ( OM0->dimension() ), xt ( IM0->dimension() ), ut ( usize) {} 
    149                 ~stateDS() {delete IM; delete OM;} 
    150                 virtual void step() { 
    151                         xt=IM->samplecond(concat ( xt,ut )); 
    152                         dt=OM->samplecond(concat ( xt,ut )); 
    153                 }; 
    154                  
    155                 virtual void log_add ( logger &L ) { 
    156                         DS::log_add ( L ); 
    157                         L_xt=L.add(IM->_rv(),"true"); 
    158                 } 
    159                 virtual void logit ( logger &L ) { 
    160                         DS::logit ( L ); 
    161                         L.logit ( L_xt,xt); 
    162                 } 
    163  
    164         }; 
     179}; 
    165180 
    166181}; //namespace 
  • bdm/stat/libEF.h

    r270 r271  
    297297                if ( evalll ) {last_lognc=est.lognc();} 
    298298        } 
    299         const epdf& _epdf() const {return est;}; 
     299        const epdf& posterior() const {return est;}; 
    300300        const eDirich* _e() const {return &est;}; 
    301301        void set_parameters ( const vec &beta0 ) { 
  • bdm/stat/loggers.h

    r270 r271  
    4848        void logit ( int id, const vec &v ) { 
    4949                it_assert_debug(id<vectors.length(),"Logger was not initialized, run init()."); 
    50                 vectors ( id ).set_row ( ind,v );} 
     50                if(id>0){ vectors ( id ).set_row ( ind,v );} 
     51        } 
    5152        //! Save values into an itfile named after \c fname. 
    5253        void itsave(const char* fname); 
  • bdm/stat/loggers_ui.h

    r263 r271  
     1 
    12/*! 
    23  \file 
     
    1920using namespace bdm; 
    2021 
     22/*! \brief UI for class RV (description of data vectors) 
     23 
     24\code 
     25rv = { 
     26        type = "rv"; //identifier of the description 
     27        // UNIQUE IDENTIFIER same names = same variable 
     28        names = ["a", "b", "c", ...];   // which will be used e.g. in loggers  
     29 
     30        //optional arguments 
     31        sizes = [1, 2, 3, ...];         // (optional) default = ones() 
     32        times = [-1, -2, 0, ...];       // time shifts with respect to current time (optional) default = zeros() 
     33} 
     34\endcode 
     35*/ 
     36class UIrv: public UIbuilder{ 
     37        public: 
     38                UIrv():UIbuilder("rv"){}; 
     39                bdmroot* build(Setting &S) const{ 
     40                        Array<string> A=get_as(S["names"]); 
     41                        ivec szs; 
     42                        ivec tms; 
     43                        if (S.exists("sizes")){  
     44                                szs=getivec(S["sizes"]); 
     45                        } else { 
     46                                szs = ones_i(A.length()); 
     47                        } 
     48                        if (S.exists("times")){  
     49                                tms=getivec(S["times"]); 
     50                        } else { 
     51                                tms = zeros_i(A.length()); 
     52                        } 
     53                        RV *tmp = new RV(A,szs,tms); 
     54                        return tmp; 
     55                }; 
     56}; 
     57UIREGISTER(UIrv); 
     58 
     59/*! \brief UI for dirfilelog (Kst file format)  
     60\code 
     61logger = { 
     62        type = "dirfilelog"; 
     63        dirmane = "directory_for_files"; // resulting files will be stored there 
     64        maxlen = 100;                    // size of memory buffer, when full results are written to disk 
     65} 
     66\endcode 
     67*/ 
    2168class UIdirfilelog : public UIbuilder { 
    2269        public: 
  • bdm/uibuilder.h

    r270 r271  
    1414#define UIREGISTER(UI) UI* UI##_global_instance = new UI(); 
    1515 
     16#define UICATCH         catch ( SettingTypeException e ) {it_error ( "Setting " +string ( e.getPath() ) +" is of incorrect Type" );}    catch ( SettingNotFoundException e ) {it_error ( "Setting " + string ( e.getPath() ) +" was not found" );} 
     17 
    1618////////// GLOBAL VAriables 
    1719 
     
    2527                UIFile ( const char * fname ) :Config() { 
    2628                        try{Config::readFile ( fname );} 
     29                        catch ( FileIOException f ) {it_error ( "File " + string ( fname ) + " not found" );} 
    2730                        catch ( ParseException& P ) { 
    2831                                char msg[200]; 
     
    3033                                it_error ( msg ); 
    3134                        } 
    32                         catch ( FileIOException f ) {it_error ( "File " + string ( fname ) + " not found" );} 
    3335                } 
    3436        }; 
     
    9597                        return tmp; 
    9698                }; 
    97         public: 
     99                const Array<string> get_as ( Setting& S ) const { 
     100                        CHECK_UITYPE ( S,TypeArray ); 
     101                        Array<string> tmp; 
     102                        tmp.set_size ( S.getLength() ); 
     103                        for ( int i=0;i<S.getLength();i++ ) {tmp(i)=(const char*)S[i];} 
     104                        return tmp; 
     105                }; 
     106                public: 
    98107                //!Constructor needs to be run only once macro UIREGISTER 
    99108                UIbuilder ( const string &typ ) {__uimap__.insert ( make_pair ( typ,this ) );} 
     
    135144                        ret = dynamic_cast<T*> ( iter->second->build ( S ) ); 
    136145                } 
    137                 catch ( SettingTypeException e ) { 
    138                         UI_DBG(S,"");  
    139                         it_error ( "Setting " +string ( e.getPath() ) +" is of incorrect Type" );} 
    140                 catch ( SettingNotFoundException e ) { 
    141                         UI_DBG(S,""); 
    142                         it_error ( "Setting " + string ( e.getPath() ) +" was not found" );} 
     146                UICATCH; 
    143147        }; 
    144148