bdmbase.h

Go to the documentation of this file.

#ifndef BDMBASE_H
#define BDMBASE_H

#include <map>

#include "../itpp_ext.h"
#include "../bdmroot.h"
#include "../shared_ptr.h"
#include "user_info.h"

using namespace libconfig;
using namespace itpp;
using namespace std;

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

class RV : public root {
private:
        typedef std::map<string, int> str2int_map;

        static ivec SIZES;
        static Array<string> NAMES;

        const static int BUFFER_STEP;

        static str2int_map MAP;

public:

protected:
        int dsize;
        int len;
        ivec ids;
        ivec times;

private:
        enum enum_dummy {dummy};

        void init ( const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times );
        int init ( const string &name, int size );
        explicit RV ( const ivec &ids0, enum_dummy dum ) : dsize ( 0 ), len ( ids0.length() ), ids ( ids0 ), times ( zeros_i ( ids0.length() ) ) {
                dsize = countsize();
        }
public:

        RV ( const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times ) {
                init ( in_names, in_sizes, in_times );
        }

        RV ( const Array<std::string> &in_names, const ivec &in_sizes ) {
                init ( in_names, in_sizes, zeros_i ( in_names.length() ) );
        }

        RV ( const 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 );

        // compiler-generated copy constructor is used


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

        string to_string() {ostringstream o; o << this; return o.str();}
        
        int _dsize() const {
                return dsize;
        }

        const ivec& _ids() const {
                return ids;
        }

        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 SIZES ( ids ( at ) );
        }
        int time ( int at ) const {
                return times ( at );
        }
        std::string name ( int at ) const {
                return NAMES ( ids ( at ) );
        }
        std::string scalarname ( int scalat ) const {
                bdm_assert(scalat<dsize,"Wrong input index");
                int id=0;
                int scalid=0;
                while (scalid+SIZES(ids(id))<=scalat)  { scalid+=SIZES(ids(id)); id++;};
                //now id is the id of variable of interest
                if (size(id)==1)
                        return  NAMES ( ids ( id ) );
                else
                        return  NAMES ( ids ( id ) )+ "_" + num2str(scalat-scalid);

        }
        void set_time ( int at, int time0 ) {
                times ( at ) = time0;
        }


        ivec findself ( const RV &rv2 ) const;
        ivec findself_ids ( 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;

        void t_plus ( int delta );

        RV remove_time() const {
                return RV ( unique ( ids ), dummy );
        }
        RV copy_t(int dt) const {
                RV tmp=*this;
                tmp.t_plus(dt);
                return tmp;
        }
        RV expand_delayes() const {
                RV rvt = this->remove_time(); //rv at t=0
                RV tmp = rvt;
                int td = mint();
                for ( int i = -1; i >= td; i-- ) {
                        rvt.t_plus ( -1 );
                        tmp.add ( rvt ); //shift u1
                }
                return tmp;
        }


        str tostr() const;
        ivec dataind ( const RV &crv ) const;
        ivec dataind_part ( const RV &crv ) const;
        void dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const;
        int mint() const {
                        return times.length()>0 ? min (times) : 0;
        }

        void from_setting ( const Setting &set );

        // TODO dodelat void to_setting( Setting &set ) const;

        static void clear_all();
        string show_all();

};
UIREGISTER ( RV );
SHAREDPTR ( RV );

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

class fnc : public root {
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 root {
protected:
        int dim;
        RV rv;

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


        virtual vec sample() const {
                bdm_error ( "not implemented" );
                return vec();
        }

        virtual mat sample_m ( int N ) const;

        virtual double evallog ( const vec &val ) const {
                bdm_error ( "not implemented" );
                return 0.0;
        }

        virtual vec evallog_m ( const mat &Val ) const;

        virtual vec evallog_m ( const Array<vec> &Avec ) const;

        virtual shared_ptr<mpdf> condition ( const RV &rv ) const;

        virtual shared_ptr<epdf> marginal ( const RV &rv ) const;

        virtual vec mean() const {
                bdm_error ( "not implemneted" );
                return vec();
        }

        virtual vec variance() const {
                bdm_error ( "not implemneted" );
                return vec();
        }

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

        bool isnamed() const {
                bool b = ( dim == rv._dsize() );
                return b;
        }

        const RV& _rv() const {
                //bdm_assert_debug ( isnamed(), "" );
                return rv;
        }


        int dimension() const {
                return dim;
        }
        void from_setting ( const Setting &set ) {
                shared_ptr<RV> r = UI::build<RV> ( set, "rv", UI::optional );
                if ( r ) {
                        set_rv ( *r );
                }
        }

};
SHAREDPTR ( epdf );

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

protected:
        void set_ep ( epdf &iepdf ) {
                ep = &iepdf;
        }
        void set_ep ( epdf *iepdfp ) {
                ep = iepdfp;
        }

public:

        mpdf() : dimc ( 0 ), rvc(), ep ( NULL ) { }

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


        virtual vec samplecond ( const vec &cond ) {
                bdm_error ( "Not implemented" );
                return vec();
        }

        virtual mat samplecond_m ( const vec &cond, int N );

        virtual double evallogcond ( const vec &dt, const vec &cond ) {
                bdm_error ( "Not implemented" );
                return 0.0;
        }

        virtual vec evallogcond_m ( const mat &Dt, const vec &cond ) {
                vec v ( Dt.cols() );
                for ( int i = 0; i < Dt.cols(); i++ ) {
                        v ( i ) = evallogcond ( Dt.get_col ( i ), cond );
                }
                return v;
        }

        virtual vec evallogcond_m ( const Array<vec> &Dt, const vec &cond ) {
                bdm_error ( "Not implemented" );
                return vec();
        }


        const RV& _rv() const {
                return ep->_rv();
        }
        const RV& _rvc() const {
                return rvc;
        }

        int dimension() const {
                return ep->dimension();
        }
        int dimensionc() {
                return dimc;
        }

        void from_setting ( const Setting &set );

        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() );
        }
};
SHAREDPTR ( mpdf );

template <class EPDF>
class mpdf_internal: public mpdf {
protected :
        EPDF iepdf;
public:
        mpdf_internal() : mpdf(), iepdf() {
                set_ep ( iepdf );
        }

        virtual void condition ( const vec &cond ) {
                bdm_error ( "Not implemented" );
        }

        EPDF& e() {
                return iepdf;
        }

        vec samplecond ( const vec &cond );
        double evallogcond ( const vec &val, const vec &cond );
        virtual vec evallogcond_m ( const mat &Dt, const vec &cond );
        virtual vec evallogcond_m ( const Array<vec> &Dt, const vec &cond );
        virtual mat samplecond_m ( const vec &cond, int N );
};

class datalink {
protected:
        int downsize;

        int upsize;

        ivec v2v_up;

public:
        datalink() : downsize ( 0 ), upsize ( 0 ) { }

        datalink ( const RV &rv, const RV &rv_up ) {
                set_connection ( rv, rv_up );
        }

        virtual void set_connection ( const RV &rv, const RV &rv_up );

        virtual void set_connection ( int ds, int us, const ivec &upind );

        vec pushdown ( const vec &val_up ) {
                vec tmp ( downsize );
                filldown ( val_up, tmp );
                return tmp;
        }
        virtual void filldown ( const vec &val_up, vec &val_down ) {
                bdm_assert_debug ( upsize == val_up.length(), "Wrong val_up" );
                val_down = val_up ( v2v_up );
        }
        virtual void pushup ( vec &val_up, const vec &val ) {
                bdm_assert_debug ( downsize == val.length(), "Wrong val" );
                bdm_assert_debug ( upsize == val_up.length(), "Wrong val_up" );
                set_subvector ( val_up, v2v_up, val );
        }
        int _upsize() {
                return upsize;
        }
        int _downsize() {
                return downsize;
        }
};

class datalink_part : public datalink {
protected:
        ivec v2v_down;
public:
        void set_connection ( const RV &rv, const RV &rv_up );
        void filldown ( const vec &val_up, vec &val_down ) {
                set_subvector ( val_down, v2v_down, val_up ( v2v_up ) );
        }
};

class datalink_buffered: public datalink_part {
protected:
        vec history;
        RV Hrv;
        ivec h2v_down;
        ivec h2v_hist;
        ivec v2h_up;
public:

        datalink_buffered() : datalink_part(), history ( 0 ), h2v_down ( 0 ), h2v_hist ( 0 ) {};
        void step ( const vec &val_up ) {
                if ( v2h_up.length() > 0 ) {
                        history.shift_right ( 0, v2h_up.length() );
                  history.set_subvector ( 0, val_up(v2h_up) ); 
                }
        }
        vec pushdown ( const vec &val_up ) {
                vec tmp ( downsize );
                filldown ( val_up, tmp );
                return tmp;
        }

        void filldown ( const vec &val_up, vec &val_down ) {
                bdm_assert_debug ( val_down.length() >= downsize, "short val_down" );

                set_subvector ( val_down, v2v_down, val_up ( v2v_up ) ); // copy direct values
                set_subvector ( val_down, h2v_down, history ( h2v_hist ) ); // copy delayed values
        }

        void set_connection ( const RV &rv, const RV &rv_up ) {
                // create link between up and down
                datalink_part::set_connection ( rv, rv_up );

                // create rvs of history
                // we can store only what we get in rv_up - everything else is removed
                ivec valid_ids = rv.findself_ids ( rv_up );
                RV rv_hist = rv.subselect ( find ( valid_ids >= 0 ) ); // select only rvs that are in rv_up
                RV rv_hist0 =rv_hist.remove_time();  // these RVs will form history at time =0
                // now we need to know what is needed from Up
                rv_hist = rv_hist.expand_delayes(); // full regressor - including time 0
                Hrv=rv_hist.subt(rv_hist0);   // remove time 0
                history = zeros ( Hrv._dsize() );

                // decide if we need to copy val to history
                if (Hrv._dsize() >0) {
                        v2h_up = rv_hist0.dataind(rv_up); // indeces of elements of rv_up to be copied
                } // else v2h_up is empty
                
                Hrv.dataind ( rv, h2v_hist, h2v_down );

                downsize = v2v_down.length() + h2v_down.length();
                upsize = v2v_up.length();
        }
        void set_history(const RV& rv1, const vec &hist0){
                bdm_assert(rv1._dsize()==hist0.length(),"hist is not compatible with given rv1");
                ivec ind_H;
                ivec ind_h0;
                Hrv.dataind(rv1, ind_H, ind_h0); // find indeces of rv in 
                set_subvector(history, ind_H, hist0(ind_h0)); // copy given hist to appropriate places
        }
};

class datalink_2to1_buffered {
protected:
        datalink_buffered dl1;
        datalink_buffered dl2;
public:
        void set_connection ( const RV &rv, const RV &rv_up1, const RV &rv_up2 ) {
                dl1.set_connection ( rv, rv_up1 );
                dl2.set_connection ( rv, rv_up2 );
        }
        void filldown ( const vec &val1, const vec &val2, vec &val_down ) {
                bdm_assert_debug ( val_down.length() >= dl1._downsize() + dl2._downsize(), "short val_down" );
                dl1.filldown ( val1, val_down );
                dl2.filldown ( val2, val_down );
        }
        void step ( const vec &dt, const vec &ut ) {
                dl1.step ( dt );
                dl2.step ( ut );
        }
};



class datalink_m2e: public datalink {
protected:
        int condsize;

        ivec v2c_up;

        ivec v2c_lo;

public:
        datalink_m2e() : condsize ( 0 ) { }

        void set_connection ( const RV &rv, const RV &rvc, const RV &rv_up );

        vec get_cond ( const vec &val_up );

        void pushup_cond ( vec &val_up, const vec &val, const vec &cond );
};

class datalink_m2m: public datalink_m2e {
protected:
        ivec c2c_up;
        ivec c2c_lo;

public:
        datalink_m2m() {};
        void set_connection ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) {
                datalink_m2e::set_connection ( rv, rvc, rv_up );
                //establish c2c connection
                rvc.dataind ( rvc_up, c2c_lo, c2c_up );
                bdm_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 root {
protected:
        Array<RV> entries;
        Array<string> names;
public:
        logger() : entries ( 0 ), names ( 0 ) {}

        virtual int add ( const RV &rv, string prefix = "" ) {
                int id;
                if ( rv._dsize() > 0 ) {
                        id = entries.length();
                        names = concat ( names, prefix ); // 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 ) {
                bdm_error ( "Not implemented" );
        };
        virtual void logit ( int id, const double &d ) {
                bdm_error ( "Not implemented" );
        };

        virtual void step() {
                bdm_error ( "Not implemneted" );
        };

        virtual void finalize() {};

        virtual void init() {};

};

class mepdf : public mpdf {
        shared_ptr<epdf> iepdf;
public:
        mepdf() { }
        mepdf ( shared_ptr<epdf> em ) {
                iepdf = em;
                set_ep ( *iepdf.get() );
                dimc = 0;
        }

        vec samplecond ( const vec &cond ) {
                return iepdf->sample();
        }
        double evallogcond ( const vec &val, const vec &cond ) {
                return iepdf->evallog ( val );
        }

        void from_setting ( const Setting &set );
};
UIREGISTER ( mepdf );
SHAREDPTR ( mepdf );

RV get_composite_rv ( const Array<shared_ptr<mpdf> > &mpdfs, bool checkoverlap = false );

class DS : public root {
protected:
        int dtsize;
        int utsize;
                int ytsize;
        RV Drv;
        RV Urv; //
                RV Yrv; //
        ivec LIDs;
public:
                DS() : Drv(), Urv(),Yrv(), LIDs(2) {};

        virtual int max_length() {return std::numeric_limits< int >::max();}
        virtual void getdata ( vec &dt ) {
                bdm_error ( "abstract class" );
        }
        virtual void getdata ( vec &dt, const ivec &indeces ) {
                bdm_error ( "abstract class" );
        }

        virtual void write (const vec &ut ) {
                bdm_error ( "abstract class" );
        }

        virtual void write (const vec &ut, const ivec &indeces ) {
                bdm_error ( "abstract class" );
        }

        virtual void step() 
        {
                bdm_error ( "abstract class" );
        }

        virtual void log_add ( logger &L ) {
                bdm_assert ( dtsize == Drv._dsize(), "invalid DS: dtsize (" + num2str ( dtsize ) + ") different from Drv " + num2str ( Drv._dsize() ) );
                bdm_assert ( utsize == Urv._dsize(), "invalid DS: utsize (" + num2str ( utsize ) + ") different from Urv " + num2str ( Urv._dsize() ) );

                LIDs(0) = L.add ( Drv, "" );
                LIDs(1) = L.add ( Urv, "" );
        }
        virtual void logit ( logger &L ) {
                vec tmp ( Drv._dsize() + Urv._dsize() );
                getdata ( tmp );
                // d is first in getdata
                L.logit ( LIDs(0), tmp.left ( Drv._dsize() ) );
                // u follows after d in getdata
                L.logit ( LIDs(1), tmp.mid ( Drv._dsize(), Urv._dsize() ) );
        }
        virtual const RV& _drv() const {
                //              return concat (Drv, Urv);// why!!!
                return Drv;// why!!!
        }
        const RV& _urv() const {
                return Urv;
        }
                const RV& _yrv() const {
                        return Yrv;
                }
                virtual void set_drv (const  RV &yrv, const RV &urv) {
                    Yrv = yrv;
                        Drv = concat(yrv,urv);
                Urv = urv;
        }
};

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

        BM() : ll ( 0 ), evalll ( true ), LIDs ( 4 ), LFlags ( 4 ) {
                LIDs = -1;/*empty IDs*/
                LFlags = 0;
                LFlags ( 0 ) = 1;  /*log only mean*/
        };
        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 {
                bdm_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 {
                bdm_error ( "Not implemented" );
                return NULL;
        };
        virtual mpdf* predictor() const {
                bdm_error ( "Not implemented" );
                return NULL;
        };


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

        virtual void condition ( const vec &val ) {
                bdm_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 void set_options ( const string &opt ) {
                LFlags ( 0 ) = 1;
                if ( opt.find ( "logbounds" ) != string::npos ) {
                        LFlags ( 1 ) = 1;
                        LFlags ( 2 ) = 1;
                }
                if ( opt.find ( "logll" ) != string::npos ) {
                        LFlags ( 3 ) = 1;
                }
        }
        ivec LIDs;

        ivec LFlags;
        virtual 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
                if ( LFlags ( 0 ) ) LIDs ( 0 ) = L.add ( r, name + "mean_" );
                if ( LFlags ( 1 ) ) LIDs ( 1 ) = L.add ( r, name + "lb_" );
                if ( LFlags ( 2 ) ) LIDs ( 2 ) = L.add ( r, name + "ub_" );
                if ( LFlags ( 3 ) ) LIDs ( 3 ) = L.add ( RV ( "ll", 1 ), name );    //TODO: "local" RV
        }
        virtual void logit ( logger &L ) {
                L.logit ( LIDs ( 0 ), posterior().mean() );
                if ( LFlags ( 1 ) || LFlags ( 2 ) ) {  //if one of them is off, its LID==-1 and will not be stored
                        vec ub, lb;
                        posterior().qbounds ( lb, ub );
                        L.logit ( LIDs ( 1 ), lb );
                        L.logit ( LIDs ( 2 ), ub );
                }
                if ( LFlags ( 3 ) ) L.logit ( LIDs ( 3 ), ll );
        }
        void from_setting ( const Setting &set ) {
                shared_ptr<RV> r = UI::build<RV> ( set, "drv", UI::optional );
                if ( r ) {
                        set_drv ( *r );
                }
                string opt;
                if ( !UI::get ( opt, set, "options", UI::optional ) ) {
                        set_options ( opt );
                }
        }

};
typedef Array<shared_ptr<epdf> > epdf_array;
typedef Array<shared_ptr<mpdf> > mpdf_array;

template<class EPDF>
vec mpdf_internal<EPDF>::samplecond ( const vec &cond ) {
        condition ( cond );
        vec temp = iepdf.sample();
        return temp;
}

template<class EPDF>
mat mpdf_internal<EPDF>::samplecond_m ( const vec &cond, int N ) {
        condition ( cond );
        mat temp ( dimension(), N );
        vec smp ( dimension() );
        for ( int i = 0; i < N; i++ ) {
                smp = iepdf.sample();
                temp.set_col ( i, smp );
        }

        return temp;
}

template<class EPDF>
double mpdf_internal<EPDF>::evallogcond ( const vec &dt, const vec &cond ) {
        double tmp;
        condition ( cond );
        tmp = iepdf.evallog ( dt );
        return tmp;
}

template<class EPDF>
vec mpdf_internal<EPDF>::evallogcond_m ( const mat &Dt, const vec &cond ) {
        condition ( cond );
        return iepdf.evallog_m ( Dt );
}

template<class EPDF>
vec mpdf_internal<EPDF>::evallogcond_m ( const Array<vec> &Dt, const vec &cond ) {
        condition ( cond );
        return iepdf.evallog_m ( Dt );
}

}; //namespace
#endif // BDMBASE_H

---