libBM.h

Go to the documentation of this file.

#ifndef BM_H
#define BM_H


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

namespace bdm {
using namespace itpp;
using namespace std;

class bdmroot {
public:
        virtual ~bdmroot() {}
};

typedef std::map<string, int> RVmap;
extern ivec RV_SIZES;
extern Array<string> RV_NAMES;

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 dsize;
        int len;
        ivec ids;
        ivec times;

private:
        void init ( Array<std::string> in_names, ivec in_sizes, ivec in_times );
        int init ( const  string &name, int size );
public:

        RV ( Array<std::string> in_names, ivec in_sizes, ivec in_times ) {init ( in_names,in_sizes,in_times );};
        RV ( Array<std::string> in_names, ivec in_sizes ) {init ( in_names,in_sizes,zeros_i ( in_names.length() ) );};
        RV ( Array<std::string> in_names ) {init ( in_names,ones_i ( in_names.length() ),zeros_i ( in_names.length() ) );}
        RV () :dsize ( 0 ),len ( 0 ),ids ( 0 ),times ( 0 ) {};
        RV ( string name, int sz, int tm=0 );


        friend std::ostream &operator<< ( std::ostream &os, const RV &rv );
        int _dsize() const {return dsize;} ;
        int countsize() const;
        ivec cumsizes() const;
        int length() const {return len;} ;
        int id ( int at ) const{return ids ( at );};
        int size ( int at ) const {return RV_SIZES ( ids ( at ) );};
        int time ( int at ) const{return times ( at );};
        std::string name ( int at ) const {return RV_NAMES ( ids ( at ) );};
        void set_time ( int at, int time0 ) {times ( at ) =time0;};

        //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 {return subselect ( ind );};
        RV operator() ( int di1, int di2 ) const {
                ivec sz=cumsizes();
                int i1=0;
                while ( sz ( i1 ) <di1 ) i1++;
                int i2=i1;
                while ( sz ( i2 ) <di2 ) i2++;
                return subselect ( linspace ( i1,i2 ) );
        };
        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 () const {return min ( times );};

};


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

extern RV RV0;


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

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

        int dimension() const{return dimy;}
};

class mpdf;


class epdf :public bdmroot {
protected:
        int dim;
        RV rv;

public:
        epdf() :dim ( 0 ),rv ( ) {};
        epdf ( const epdf &e ) :dim ( e.dim ),rv ( e.rv ) {};
        epdf ( const RV &rv0 ) {set_rv ( rv0 );};
        void set_parameters ( int dim0 ) {dim=dim0;}


        virtual vec sample () const {it_error ( "not implemneted" );return vec ( 0 );};
        virtual mat sample_m ( int N ) const;
        virtual double evallog ( const vec &val ) const {it_error ( "not implemneted" );return 0.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 {it_error ( "not implemneted" );return vec ( 0 );};
        virtual vec variance() const {it_error ( "not implemneted" );return vec ( 0 );};
        virtual void qbounds ( vec &lb, vec &ub, double percentage=0.95 ) const {
                vec mea=mean(); vec std=sqrt(variance()); 
                lb = mea-2*std; ub=mea+2*std;
        };


        void set_rv ( const RV &rv0 ) {rv = rv0; }//it_assert_debug(isnamed(),""); };
        bool isnamed() const {bool b= ( dim==rv._dsize() );return b;}
        const RV& _rv() const {it_assert_debug ( isnamed(),"" ); return rv;}


        int dimension() const {return dim;}

};


//TODO Samplecond can be generalized

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

        mpdf ( ) :dimc ( 0 ),rvc ( ) {};
        mpdf ( const mpdf &m ) :dimc ( m.dimc ),rvc ( m.rvc ) {};


        virtual vec samplecond ( const vec &cond ) {
                this->condition ( cond );
                vec temp= ep->sample();
                return temp;
        };
        virtual mat samplecond_m ( const vec &cond, int N ) {
                this->condition ( cond );
                mat temp ( ep->dimension(),N ); vec smp ( ep->dimension() );
                for ( int i=0;i<N;i++ ) {smp=ep->sample() ;temp.set_col ( i, 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 );};


        RV _rv() {return ep->_rv();}
        RV _rvc() {it_assert_debug ( isnamed(),"" ); return rvc;}
        int dimension() {return ep->dimension();}
        int dimensionc() {return dimc;}
        epdf& _epdf() {return *ep;}
        epdf* _e() {return ep;}

        void set_rvc ( const RV &rvc0 ) {rvc=rvc0;}
        void set_rv ( const RV &rv0 ) {ep->set_rv ( rv0 );}
        bool isnamed() {return ( ep->isnamed() ) && ( dimc==rvc._dsize() );}
};

class datalink {
protected:
        int downsize;
        int upsize;
        ivec v2v_up;
public:
        datalink () {};
        datalink ( const RV &rv, const RV &rv_up ) {set_connection ( rv,rv_up );};
        void set_connection ( const RV &rv, const RV &rv_up ) {
                downsize = rv._dsize();
                upsize = rv_up._dsize();
                v2v_up= ( rv.dataind ( rv_up ) );

                it_assert_debug ( v2v_up.length() ==downsize,"rv is not fully in rv_up" );
        }
        vec pushdown ( const vec &val_up ) {
                it_assert_debug ( upsize==val_up.length(),"Wrong val_up" );
                return get_vec ( val_up,v2v_up );
        }
        void pushup ( vec &val_up, const vec &val ) {
                it_assert_debug ( downsize==val.length(),"Wrong val" );
                it_assert_debug ( upsize==val_up.length(),"Wrong val_up" );
                set_subvector ( val_up, v2v_up, val );
        }
};

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

public:
        datalink_m2e ( const RV &rv,  const RV &rvc, const RV &rv_up ) :
                        datalink ( rv,rv_up ), condsize ( rvc._dsize() ) {
                //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 pushup_cond ( vec &val_up, const vec &val, const vec &cond ) {
                it_assert_debug ( downsize==val.length(),"Wrong val" );
                it_assert_debug ( upsize==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 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;
                if ( rv._dsize() >0 ) {
                        id=entries.length();
                        names=concat ( names, name ); // diff
                        entries.set_length ( id+1,true );
                        entries ( id ) = rv;
                }
                else { id =-1;}
                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() {};

};

class mepdf : public mpdf {
public:
        mepdf ( const epdf* em ) :mpdf ( ) {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:
        int dtsize;
        int utsize;
        RV Drv;
        RV Urv; //
        int L_dt, L_ut;
public:
        DS() :Drv ( ),Urv ( ) {};
        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 ) {
                it_assert_debug ( dtsize==Drv._dsize(),"" );
                it_assert_debug ( utsize==Urv._dsize(),"" );

                L_dt=L.add ( Drv,"" );
                L_ut=L.add ( Urv,"" );
        }
        virtual void logit ( logger &L ) {
                vec tmp ( Drv._dsize() +Urv._dsize() );
                getdata ( tmp );
                // d is first in getdata
                L.logit ( L_dt,tmp.left ( Drv._dsize() ) );
                // u follows after d in getdata
                L.logit ( L_ut,tmp.mid ( Drv._dsize(), Urv._dsize() ) );
        }
        virtual RV _drv() const {return concat ( Drv,Urv );}
        const RV& _urv() const {return Urv;}
};

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

        BM () :ll ( 0 ),evalll ( true ), LIDs ( 3 ), opt_L_bounds ( false ) {};
        BM ( const BM &B ) :  drv ( B.drv ), ll ( B.ll ), evalll ( B.evalll ) {}
        virtual BM* _copy_ () const {return NULL;};


        virtual void bayes ( const vec &dt ) = 0;
        virtual void bayesB ( const mat &Dt );
        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* epredictor ( ) const {it_error ( "Not implemented" );return NULL;};
        virtual mpdf* predictor ( ) const {it_error ( "Not implemented" );return NULL;};


        RV rvc;
        const RV& _rvc() const {return rvc;}

        virtual void condition ( const vec &val ) {it_error ( "Not implemented!" );};




        const RV& _drv() const {return drv;}
        void set_drv ( const RV &rv ) {drv=rv;}
        void set_rv ( const RV &rv ) {const_cast<epdf&> ( posterior() ).set_rv ( rv );}
        double _ll() const {return ll;}
        void set_evalll ( bool evl0 ) {evalll=evl0;}
        virtual const epdf& posterior() const =0;
        virtual const epdf* _e() const =0;


        void set_options ( const string &opt ) {
                opt_L_bounds= ( opt.find ( "logbounds" ) !=string::npos );
        }
        ivec LIDs;

        bool opt_L_bounds;
        void log_add ( logger *L, const string &name="" ) {
                // internal
                RV r;
                if ( posterior().isnamed() ) {r=posterior()._rv();}
                else{r=RV ( "est", posterior().dimension() );};

                // Add mean value
                LIDs ( 0 ) =L->add ( r,name );
                if ( opt_L_bounds ) {
                        LIDs ( 1 ) =L->add ( r,name+"_lb" );
                        LIDs ( 2 ) =L->add ( r,name+"_ub" );
                }
        }
        void logit ( logger *L ) {
                L->logit ( LIDs ( 0 ), posterior().mean() );
                if ( opt_L_bounds ) {
                        vec ub,lb;
                        posterior().qbounds(lb,ub);
                        L->logit ( LIDs ( 1 ), lb );
                        L->logit ( LIDs ( 2 ), ub );
                }
        }
};

}; //namespace
#endif // BM_H

---