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;
        int length() const {return len;} ;
        int id ( int at ) const{return ids ( at );};
        int size ( int at ) const {return RV_SIZES ( at );};
        int time ( int at ) const{return times ( at );};
        std::string name ( int at ) const {return RV_NAMES ( 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 );};
        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 _dimy() 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) {};
        
        
        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;
        
        
        void set_rv ( const RV &rv0 ) {rv = rv0; it_assert_debug(isnamed(),""); };
        bool isnamed() const {return ( dim==rv._dsize() );}
        const RV& _rv() const {it_assert_debug ( isnamed(),"" ); return rv;}
        
        
        bool 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 ( 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=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() {};

};

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 ( false) {};
        BM ( const BM &B ) :  drv(B.drv), ll ( B.ll ), evalll ( B.evalll ) {}
        virtual BM* _copy_ () {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;};
        
        
        const RV& _drv() const {return drv;}
        void set_drv ( const RV &rv ) {drv=rv;}
        double _ll() const {return ll;}
        void set_evalll ( bool evl0 ) {evalll=evl0;}
        virtual const epdf& _epdf() const =0;
        virtual const epdf* _e() const =0;

};

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

}; //namespace
#endif // BM_H

---