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

Go to the documentation of this file.

#ifndef BM_H
#define BM_H

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

using namespace itpp;

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

        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 {
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 {
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 ~epdf() {};
        const RV& _rv() const {return rv;}
        void _renewrv ( const RV &in_rv ) {rv=in_rv;}
};


//TODO Samplecond can be generalized

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

//      virtual fnc moment ( const int order = 1 );
        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;}
};

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 {
protected:
        RV Drv;
        RV Urv; //
public:
        void getdata ( vec &dt );
        void getdata ( vec &dt, ivec &indeces );
        void write ( vec &ut );
        void write ( vec &ut, ivec &indeces );
        void linkrvs ( RV &drv, RV &urv );

        void step();

};

class BM {
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 {
protected:
        RV rvc;
public:
        virtual void condition ( const vec &val ) =0;
        BMcond ( RV &rv0 ) :rvc ( rv0 ) {};
        virtual ~BMcond() {};
        const RV& _rvc() const {return rvc;}
};

#endif // BM_H

---