- Main Page
- Related Pages
- Modules
- Classes
- Files

# work/git/mixpp/bdm/stat/libBM.h

Go to the documentation of this file.

    00001
    00017 #ifndef BM_H
    00018 #define BM_H
    00019
    00020 #include <itpp/itbase.h>
    00021 #include "../itpp_ext.h"
    00022 //#include <std>
    00023
    00024 using namespace itpp;
    00025
    00027 class str {
    00028 public:
    00030         ivec ids;
    00032         ivec times;
    00034         str ( ivec ids0, ivec times0 ) :ids ( ids0 ),times ( times0 ) {
    00035                 it_assert_debug ( times0.length() ==ids0.length(),"Incompatible input" );
    00036         };
    00037 };
    00038
    00045 class RV {
    00046 protected:
    00048         int tsize;
    00050         int len;
    00052         ivec ids;
    00054         ivec sizes;
    00056         ivec times;
    00058         Array<std::string> names;
    00059
    00060 private:
    00062         void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times );
    00063 public:
    00065         RV ( Array<std::string> in_names, ivec in_sizes, ivec in_times );
    00067         RV ( Array<std::string> in_names, ivec in_sizes );
    00069         RV ( Array<std::string> in_names );
    00071         RV ();
    00072
    00074         friend std::ostream &operator<< ( std::ostream &os, const RV &rv );
    00075
    00077         int count() const {return tsize;} ;
    00079         int length() const {return len;} ;
    00080
    00081         //TODO why not inline and later??
    00082
    00084         ivec findself ( const RV &rv2 ) const;
    00086         bool equal ( const RV &rv2 ) const;
    00088         bool add ( const RV &rv2 );
    00090         RV subt ( const RV &rv2 ) const;
    00092         RV subselect ( const ivec &ind ) const;
    00094         RV operator() ( const ivec &ind ) const;
    00096         void t ( int delta );
    00098         str tostr() const;
    00101         ivec dataind ( const RV &crv ) const;
    00104         void dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const;
    00105
    00107         Array<std::string>& _names() {return names;};
    00108
    00110         int id ( int at ) {return ids ( at );};
    00112         int size ( int at ) {return sizes ( at );};
    00114         int time ( int at ) {return times ( at );};
    00116         std::string name ( int at ) {return names ( at );};
    00117
    00119         void set_id ( int at, int id0 ) {ids ( at ) =id0;};
    00121         void set_size ( int at, int size0 ) {sizes ( at ) =size0; tsize=sum ( sizes );};
    00123         void set_time ( int at, int time0 ) {times ( at ) =time0;};
    00124
    00126         void newids();
    00127 };
    00128
    00130 RV concat ( const RV &rv1, const RV &rv2 );
    00131
    00133 extern RV RV0;
    00134
    00136
    00137 class fnc {
    00138 protected:
    00140         int dimy;
    00141 public:
    00143         fnc ( int dy ) :dimy ( dy ) {};
    00145         virtual vec eval ( const vec &cond ) {
    00146                 return vec ( 0 );
    00147         };
    00148
    00150         int _dimy() const{return dimy;}
    00151
    00153         virtual ~fnc() {};
    00154 };
    00155
    00156 class mpdf;
    00157
    00159
    00160 class epdf {
    00161 protected:
    00163         RV rv;
    00164 public:
    00166         epdf() :rv ( ) {};
    00167
    00169         epdf ( const RV &rv0 ) :rv ( rv0 ) {};
    00170
    00171 //      //! Returns the required moment of the epdf
    00172 //      virtual vec moment ( const int order = 1 );
    00173
    00175         virtual vec sample () const =0;
    00177         virtual mat sample_m ( int N ) const;
    00178
    00180         virtual double evallog ( const vec &val ) const =0;
    00181
    00183         virtual vec evallog_m ( const mat &Val ) const {
    00184                 vec x ( Val.cols() );
    00185                 for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog ( Val.get_col ( i ) ) ;}
    00186                 return x;
    00187         }
    00189         virtual mpdf* condition ( const RV &rv ) const  {it_warning ( "Not implemented" ); return NULL;}
    00191         virtual epdf* marginal ( const RV &rv ) const {it_warning ( "Not implemented" ); return NULL;}
    00192
    00194         virtual vec mean() const =0;
    00195
    00197         virtual ~epdf() {};
    00199         const RV& _rv() const {return rv;}
    00201         void _renewrv ( const RV &in_rv ) {rv=in_rv;}
    00203 };
    00204
    00205
    00207 //TODO Samplecond can be generalized
    00208
    00209 class mpdf {
    00210 protected:
    00212         RV rv;
    00214         RV rvc;
    00216         epdf* ep;
    00217 public:
    00218
    00220 //      virtual fnc moment ( const int order = 1 );
    00222         virtual vec samplecond ( const vec &cond, double &ll ) {
    00223                 this->condition ( cond );
    00224                 vec temp= ep->sample();
    00225                 ll=ep->evallog ( temp );return temp;
    00226         };
    00228         virtual mat samplecond_m ( const vec &cond, vec &ll, int N ) {
    00229                 this->condition ( cond );
    00230                 mat temp ( rv.count(),N ); vec smp ( rv.count() );
    00231                 for ( int i=0;i<N;i++ ) {smp=ep->sample() ;temp.set_col ( i, smp );ll ( i ) =ep->evallog ( smp );}
    00232                 return temp;
    00233         };
    00235         virtual void condition ( const vec &cond ) {it_error ( "Not implemented" );};
    00236
    00238         virtual double evallogcond ( const vec &dt, const vec &cond ) {double tmp; this->condition ( cond );tmp = ep->evallog ( dt );           it_assert_debug(std::isfinite(tmp),"Infinite value"); return tmp;
    00239         };
    00240
    00242         virtual vec evallogcond_m ( const mat &Dt, const vec &cond ) {this->condition ( cond );return ep->evallog_m ( Dt );};
    00243
    00245         virtual ~mpdf() {};
    00246
    00248         mpdf ( const RV &rv0, const RV &rvc0 ) :rv ( rv0 ),rvc ( rvc0 ) {};
    00250         RV _rvc() const {return rvc;}
    00252         RV _rv() const {return rv;}
    00254         epdf& _epdf() {return *ep;}
    00255 };
    00256
    00259 class datalink_e2e {
    00260 protected:
    00262         int valsize;
    00264         int valupsize;
    00266         ivec v2v_up;
    00267 public:
    00269         datalink_e2e ( const RV &rv, const RV &rv_up ) :
    00270                         valsize ( rv.count() ), valupsize ( rv_up.count() ), v2v_up ( rv.dataind ( rv_up ) )  {
    00271                 it_assert_debug ( v2v_up.length() ==valsize,"rv is not fully in rv_up" );
    00272         }
    00274         vec get_val ( const vec &val_up ) {it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" ); return get_vec ( val_up,v2v_up );}
    00276         void fill_val ( vec &val_up, const vec &val ) {
    00277                 it_assert_debug ( valsize==val.length(),"Wrong val" );
    00278                 it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" );
    00279                 set_subvector ( val_up, v2v_up, val );
    00280         }
    00281 };
    00282
    00284 class datalink_m2e: public datalink_e2e {
    00285 protected:
    00287         int condsize;
    00289         ivec v2c_up;
    00291         ivec v2c_lo;
    00292
    00293 public:
    00295         datalink_m2e ( const RV &rv,  const RV &rvc, const RV &rv_up ) :
    00296                         datalink_e2e ( rv,rv_up ), condsize ( rvc.count() ) {
    00297                 //establish v2c connection
    00298                 rvc.dataind ( rv_up, v2c_lo, v2c_up );
    00299         }
    00301         vec get_cond ( const vec &val_up ) {
    00302                 vec tmp ( condsize );
    00303                 set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) );
    00304                 return tmp;
    00305         }
    00306         void fill_val_cond ( vec &val_up, const vec &val, const vec &cond ) {
    00307                 it_assert_debug ( valsize==val.length(),"Wrong val" );
    00308                 it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" );
    00309                 set_subvector ( val_up, v2v_up, val );
    00310                 set_subvector ( val_up, v2c_up, cond );
    00311         }
    00312 };
    00315 class datalink_m2m: public datalink_m2e {
    00316 protected:
    00318         ivec c2c_up;
    00320         ivec c2c_lo;
    00321 public:
    00323         datalink_m2m ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) :
    00324                         datalink_m2e ( rv, rvc, rv_up) {
    00325                 //establish c2c connection
    00326                 rvc.dataind ( rvc_up, c2c_lo, c2c_up );
    00327                 it_assert_debug(c2c_lo.length()+v2c_lo.length()==condsize, "cond is not fully given");
    00328         }
    00330         vec get_cond ( const vec &val_up, const vec &cond_up ) {
    00331                 vec tmp ( condsize );
    00332                 set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) );
    00333                 set_subvector ( tmp,c2c_lo,cond_up ( c2c_up ) );
    00334                 return tmp;
    00335         }
    00337
    00338 };
    00339
    00343 class mepdf : public mpdf {
    00344 public:
    00346         mepdf (const epdf* em ) :mpdf ( em->_rv(),RV() ) {ep=const_cast<epdf*>(em);};
    00347         void condition ( const vec &cond ) {}
    00348 };
    00349
    00352 class compositepdf {
    00353 protected:
    00355         int n;
    00357         Array<mpdf*> mpdfs;
    00358 public:
    00359         compositepdf ( Array<mpdf*> A0 ) : n ( A0.length() ), mpdfs ( A0 ) {};
    00361         RV getrv ( bool checkoverlap=false );
    00363         void setrvc ( const RV &rv, RV &rvc );
    00364 };
    00365
    00373 class DS {
    00374 protected:
    00376         RV Drv;
    00378         RV Urv; //
    00379 public:
    00381         void getdata ( vec &dt );
    00383         void getdata ( vec &dt, ivec &indeces );
    00385         void write ( vec &ut );
    00387         void write ( vec &ut, ivec &indeces );
    00393         void linkrvs ( RV &drv, RV &urv );
    00394
    00396         void step();
    00397
    00398 };
    00399
    00404 class BM {
    00405 protected:
    00407         RV rv;
    00409         double ll;
    00411         bool evalll;
    00412 public:
    00413
    00415         BM ( const RV &rv0, double ll0=0,bool evalll0=true ) :rv ( rv0 ), ll ( ll0 ),evalll ( evalll0 ) {//Fixme: test rv
    00416         };
    00418         BM ( const BM &B ) : rv ( B.rv ), ll ( B.ll ), evalll ( B.evalll ) {}
    00419
    00423         virtual void bayes ( const vec &dt ) = 0;
    00425         virtual void bayesB ( const mat &Dt );
    00427         virtual const epdf& _epdf() const =0;
    00428
    00430         virtual const epdf* _e() const =0;
    00431
    00434         virtual double logpred ( const vec &dt ) const{it_error ( "Not implemented" );return 0.0;}
    00436         vec logpred_m ( const mat &dt ) const{vec tmp ( dt.cols() );for ( int i=0;i<dt.cols();i++ ) {tmp ( i ) =logpred ( dt.get_col ( i ) );}return tmp;}
    00437
    00439         virtual epdf* predictor ( const RV &rv ) const {it_error ( "Not implemented" );return NULL;};
    00440
    00442         virtual ~BM() {};
    00444         const RV& _rv() const {return rv;}
    00446         double _ll() const {return ll;}
    00448         void set_evalll ( bool evl0 ) {evalll=evl0;}
    00449
    00452         virtual BM* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
    00453 };
    00454
    00464 class BMcond {
    00465 protected:
    00467         RV rvc;
    00468 public:
    00470         virtual void condition ( const vec &val ) =0;
    00472         BMcond ( RV &rv0 ) :rvc ( rv0 ) {};
    00474         virtual ~BMcond() {};
    00476         const RV& _rvc() const {return rvc;}
    00477 };
    00478
    00480 #endif // BM_H