libBM.h

Go to the documentation of this file.

#ifndef BM_H
#define BM_H


#include "../itpp_ext.h"
//#include <std>

namespace bdm {
        using namespace itpp;
        using std::string;

        class bdmroot {
                virtual void print() {}
        };

        class str {
        public:
                ivec ids;
                ivec times;
                str ( ivec ids0, ivec times0 ) :ids ( ids0 ),times ( times0 ) {
                        it_assert_debug ( times0.length() ==ids0.length(),"Incompatible input" );
                };
        };

        class RV :public bdmroot {
        protected:
                int tsize;
                int len;
                ivec ids;
                ivec sizes;
                ivec times;
                Array<std::string> names;

        private:
                void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times );
        public:
                RV ( Array<std::string> in_names, ivec in_sizes, ivec in_times );
                RV ( Array<std::string> in_names, ivec in_sizes );
                RV ( Array<std::string> in_names );
                RV ();

                friend std::ostream &operator<< ( std::ostream &os, const RV &rv );

                int count() const {return tsize;} ;
                int length() const {return len;} ;

                //TODO why not inline and later??

                ivec findself ( const RV &rv2 ) const;
                bool equal ( const RV &rv2 ) const;
                bool add ( const RV &rv2 );
                RV subt ( const RV &rv2 ) const;
                RV subselect ( const ivec &ind ) const;
                RV operator() ( const ivec &ind ) const;
                void t ( int delta );
                str tostr() const;
                ivec dataind ( const RV &crv ) const;
                void dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const;
                int mint () {return min ( times );};

                Array<std::string>& _names() {return names;};

                int id ( int at ) {return ids ( at );};
                int size ( int at ) {return sizes ( at );};
                int time ( int at ) {return times ( at );};
                std::string name ( int at ) {return names ( at );};

                void set_id ( int at, int id0 ) {ids ( at ) =id0;};
                void set_size ( int at, int size0 ) {sizes ( at ) =size0; tsize=sum ( sizes );};
                void set_time ( int at, int time0 ) {times ( at ) =time0;};

                void newids();
        };

        RV concat ( const RV &rv1, const RV &rv2 );

        extern RV RV0;


        class fnc :public bdmroot {
        protected:
                int dimy;
        public:
                fnc ( int dy ) :dimy ( dy ) {};
                virtual vec eval ( const vec &cond ) {
                        return vec ( 0 );
                };

                virtual void condition ( const vec &val ) {};

                int _dimy() const{return dimy;}

                virtual ~fnc() {};
        };

        class mpdf;


        class epdf :public bdmroot {
        protected:
                RV rv;
        public:
                epdf() :rv ( ) {};

                epdf ( const RV &rv0 ) :rv ( rv0 ) {};

//      //! Returns the required moment of the epdf
//      virtual vec moment ( const int order = 1 );

                virtual vec sample () const =0;
                virtual mat sample_m ( int N ) const;

                virtual double evallog ( const vec &val ) const =0;

                virtual vec evallog_m ( const mat &Val ) const {
                        vec x ( Val.cols() );
                        for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog ( Val.get_col ( i ) ) ;}
                        return x;
                }
                virtual mpdf* condition ( const RV &rv ) const  {it_warning ( "Not implemented" ); return NULL;}
                virtual epdf* marginal ( const RV &rv ) const {it_warning ( "Not implemented" ); return NULL;}

                virtual vec mean() const =0;

                virtual vec variance() const = 0;

                virtual ~epdf() {};
                const RV& _rv() const {return rv;}
                void _renewrv ( const RV &in_rv ) {rv=in_rv;}
        };


//TODO Samplecond can be generalized

        class mpdf : public bdmroot {
        protected:
                RV rv;
                RV rvc;
                epdf* ep;
        public:

                virtual vec samplecond ( const vec &cond, double &ll ) {
                        this->condition ( cond );
                        vec temp= ep->sample();
                        ll=ep->evallog ( temp );return temp;
                };
                virtual mat samplecond_m ( const vec &cond, vec &ll, int N ) {
                        this->condition ( cond );
                        mat temp ( rv.count(),N ); vec smp ( rv.count() );
                        for ( int i=0;i<N;i++ ) {smp=ep->sample() ;temp.set_col ( i, smp );ll ( i ) =ep->evallog ( smp );}
                        return temp;
                };
                virtual void condition ( const vec &cond ) {it_error ( "Not implemented" );};

                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;
                };

                virtual vec evallogcond_m ( const mat &Dt, const vec &cond ) {this->condition ( cond );return ep->evallog_m ( Dt );};

                virtual ~mpdf() {};

                mpdf ( const RV &rv0, const RV &rvc0 ) :rv ( rv0 ),rvc ( rvc0 ) {};
                RV _rvc() const {return rvc;}
                RV _rv() const {return rv;}
                epdf& _epdf() {return *ep;}
                epdf* _e() {return ep;}
        };

        class logger : public bdmroot {
        protected:
                Array<RV> entries;
                Array<string> names;
        public:
                logger ( ) : entries ( 0 ),names ( 0 ) {}

                virtual int add ( const RV &rv, string name="" ) {
                        int id=entries.length();
                        names=concat ( names, name ); // diff
                        entries.set_length ( id+1,true );
                        entries ( id ) = rv;
                        return id; // identifier of the last entry
                }

                virtual void logit ( int id, const vec &v ) =0;

                virtual void step() =0;

                virtual void finalize() {};

                virtual void init() {};

                virtual ~logger() {};
        };

        class datalink_e2e {
        protected:
                int valsize;
                int valupsize;
                ivec v2v_up;
        public:
                datalink_e2e ( const RV &rv, const RV &rv_up ) :
                                valsize ( rv.count() ), valupsize ( rv_up.count() ), v2v_up ( rv.dataind ( rv_up ) )  {
                        it_assert_debug ( v2v_up.length() ==valsize,"rv is not fully in rv_up" );
                }
                vec get_val ( const vec &val_up ) {it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" ); return get_vec ( val_up,v2v_up );}
                void fill_val ( vec &val_up, const vec &val ) {
                        it_assert_debug ( valsize==val.length(),"Wrong val" );
                        it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" );
                        set_subvector ( val_up, v2v_up, val );
                }
        };

        class datalink_m2e: public datalink_e2e {
        protected:
                int condsize;
                ivec v2c_up;
                ivec v2c_lo;

        public:
                datalink_m2e ( const RV &rv,  const RV &rvc, const RV &rv_up ) :
                                datalink_e2e ( rv,rv_up ), condsize ( rvc.count() ) {
                        //establish v2c connection
                        rvc.dataind ( rv_up, v2c_lo, v2c_up );
                }
                vec get_cond ( const vec &val_up ) {
                        vec tmp ( condsize );
                        set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) );
                        return tmp;
                }
                void fill_val_cond ( vec &val_up, const vec &val, const vec &cond ) {
                        it_assert_debug ( valsize==val.length(),"Wrong val" );
                        it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" );
                        set_subvector ( val_up, v2v_up, val );
                        set_subvector ( val_up, v2c_up, cond );
                }
        };
        class datalink_m2m: public datalink_m2e {
        protected:
                ivec c2c_up;
                ivec c2c_lo;
        public:
                datalink_m2m ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) :
                                datalink_m2e ( rv, rvc, rv_up ) {
                        //establish c2c connection
                        rvc.dataind ( rvc_up, c2c_lo, c2c_up );
                        it_assert_debug ( c2c_lo.length() +v2c_lo.length() ==condsize, "cond is not fully given" );
                }
                vec get_cond ( const vec &val_up, const vec &cond_up ) {
                        vec tmp ( condsize );
                        set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) );
                        set_subvector ( tmp,c2c_lo,cond_up ( c2c_up ) );
                        return tmp;
                }

        };

        class mepdf : public mpdf {
        public:
                mepdf ( const epdf* em ) :mpdf ( em->_rv(),RV() ) {ep=const_cast<epdf*> ( em );};
                void condition ( const vec &cond ) {}
        };

        class compositepdf {
        protected:
                int n;
                Array<mpdf*> mpdfs;
        public:
                compositepdf ( Array<mpdf*> A0 ) : n ( A0.length() ), mpdfs ( A0 ) {};
                RV getrv ( bool checkoverlap=false );
                void setrvc ( const RV &rv, RV &rvc );
        };

        class DS : public bdmroot {
        protected:
                RV Drv;
                RV Urv; //
                int L_dt, L_ut;
        public:
                DS() :Drv ( RV0 ),Urv ( RV0 ) {};
                DS ( const RV &Drv0, const RV &Urv0 ) :Drv ( Drv0 ),Urv ( Urv0 ) {};
                virtual void getdata ( vec &dt ) {it_error ( "abstract class" );};
                virtual void getdata ( vec &dt, const ivec &indeces ) {it_error ( "abstract class" );};
                virtual void write ( vec &ut ) {it_error ( "abstract class" );};
                virtual void write ( vec &ut, const ivec &indeces ) {it_error ( "abstract class" );};

                virtual void step() =0;

                virtual void log_add ( logger &L ) {
                        L_dt=L.add ( Drv,"" );
                        L_ut=L.add ( Urv,"" );
                }
                virtual void logit ( logger &L ) {
                        vec tmp(Drv.count()+Urv.count());
                        getdata(tmp);
                        // d is first in getdata
                        L.logit ( L_dt,tmp.left ( Drv.count() ) );
                        // u follows after d in getdata
                        L.logit ( L_ut,tmp.mid ( Drv.count(), Urv.count() ) );
                }
        };

        class BM :public bdmroot {
        protected:
                RV rv;
                double ll;
                bool evalll;
        public:

                BM ( const RV &rv0, double ll0=0,bool evalll0=true ) :rv ( rv0 ), ll ( ll0 ),evalll ( evalll0 ) {//Fixme: test rv
                };
                BM ( const BM &B ) : rv ( B.rv ), ll ( B.ll ), evalll ( B.evalll ) {}

                virtual void bayes ( const vec &dt ) = 0;
                virtual void bayesB ( const mat &Dt );
                virtual const epdf& _epdf() const =0;

                virtual const epdf* _e() const =0;

                virtual double logpred ( const vec &dt ) const{it_error ( "Not implemented" );return 0.0;}
                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;}

                virtual epdf* predictor ( const RV &rv ) const {it_error ( "Not implemented" );return NULL;};

                virtual ~BM() {};
                const RV& _rv() const {return rv;}
                double _ll() const {return ll;}
                void set_evalll ( bool evl0 ) {evalll=evl0;}

                virtual BM* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
        };

        class BMcond :public bdmroot {
        protected:
                RV rvc;
        public:
                virtual void condition ( const vec &val ) =0;
                BMcond ( RV &rv0 ) :rvc ( rv0 ) {};
                virtual ~BMcond() {};
                const RV& _rvc() const {return rvc;}
        };

}; //namespace
#endif // BM_H

---