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
{

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) {
                        bdm_assert_debug (times0.length() == ids0.length(), "Incompatible input");
                };
};

class RV : public root
{
        protected:
                int dsize;
                int len;
                ivec ids;
                ivec times;

        private:
                void init (const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times);
                int init (const string &name, int size);
        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);

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


                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;

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

                void from_setting (const Setting &set);

                // TODO dodelat void to_setting( Setting &set ) const;
                static void clear_all();
};
UIREGISTER (RV);
SHAREDPTR (RV);

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

extern RV RV0;


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


                RV _rv() const {
                        return ep->_rv();
                }
                RV _rvc() {
                        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);
                }

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

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

                vec pushdown (const vec &val_up) {
                        bdm_assert_debug (upsize == val_up.length(), "Wrong val_up");
                        return get_vec (val_up, v2v_up);
                }

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

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;
                RV Drv;
                RV Urv; //
                int L_dt, L_ut;
        public:
                DS() : Drv(), Urv() {};

                virtual void getdata (vec &dt) {
                        bdm_error ("abstract class");
                }

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

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

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

                virtual void step() = 0;

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

                        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;
                }
                virtual void set_drv (const  RV &drv, const RV &urv) {
                        Drv = drv;
                        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

---