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

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



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

        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 sampleN ( int N ) const;
        virtual double eval ( const vec &val ) const {return exp ( this->evalpdflog ( val ) );};

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

        virtual vec evalpdflog_m ( const mat &Val ) const {
                vec x ( Val.cols() );
                for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evalpdflog( Val.get_col(i) ) ;}
                return x;
        }
        mpdf* condition(const RV &rv){it_warning("Not implemented"); return NULL;}
        epdf* marginal(const RV &rv){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->evalpdflog ( temp );return temp;};
        virtual mat samplecond (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->evalpdflog ( smp );}
                return temp;
        };
        virtual void condition ( const vec &cond ) {it_error("Not implemented");};

        virtual double evalcond ( const vec &dt, const vec &cond ) {this->condition ( cond );return ep->eval ( dt );};

        virtual ~mpdf() {};

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

class mepdf : public mpdf {
public:
        mepdf (epdf &em ) :mpdf ( em._rv(),RV() ) {ep=&em;};
};

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

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

---