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

class RV : public root
{
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;
  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;};

  //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);};
  RV operator()(int di1, int di2) const {
    ivec sz = cumsizes();
    int i1 = 0;
    while (sz(i1) < di1) i1++;
    int i2 = i1;
    while (sz(i2) < di2) i2++;
    return subselect(linspace(i1, i2));
  };
  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);};

  // TODO aktualizovat dle soucasneho UI
  void from_setting(const Setting &set);

  // TODO dodelat void to_setting( Setting &set ) const;
};
UIREGISTER(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) {set_rv(rv0);};
  void set_parameters(int dim0) {dim = dim0;}


  virtual vec sample() const {it_error("not implemneted"); return vec(0);};
  virtual mat sample_m(int N) const;
  virtual double evallog(const vec &val) const {it_error("not implemneted"); return 0.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 vec evallog_m(const Array<vec> &Avec) const {
          vec x(Avec.size());
          for (int i = 0; i < Avec.size(); i++) {x(i) = evallog(Avec(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 {it_error("not implemneted"); return vec(0);};

  virtual vec variance() const {it_error("not implemneted"); return vec(0);};
  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; }   //it_assert_debug(isnamed(),""); };
  bool isnamed() const {bool b = (dim == rv._dsize()); return b;}
  const RV& _rv() const {it_assert_debug(isnamed(), ""); return rv;}


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

};


//TODO Samplecond can be generalized

class mpdf : public root
{
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);};
  virtual vec evallogcond_m(const Array<vec> &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 from_setting(const Setting &set){
          if (set.exists("rv")){
                  RV* r = UI::build<RV>(set,"rv");
                  set_rv(*r); 
                  delete r;
          }
          if (set.exists("rvc")){
                  RV* r = UI::build<RV>(set,"rvc");
                  set_rvc(*r); 
                  delete r;
          }
  }

  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() {};
  datalink(const RV &rv, const RV &rv_up) {set_connection(rv, rv_up);};
  void set_connection(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");
  }
  void set_connection(int ds, int us, const ivec &upind) {
    downsize = ds;
    upsize = us;
    v2v_up = upind;

    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() {};
  void set_connection(const RV &rv,  const RV &rvc, const RV &rv_up) {
    datalink::set_connection(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() {};
  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);
    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 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) = 0;
  virtual void logit(int id, const double &d) = 0;

  virtual void step() = 0;

  virtual void finalize() {};

  virtual void init() {};

};

class mepdf : public mpdf {
        bool owning_ep; // flag trigers deleting ep by destructor
public:
        mepdf(){};
        mepdf ( epdf* em, bool owning_ep0=false ) :mpdf ( ) {ep= em ;owning_ep=owning_ep0;dimc=0;};
        mepdf (const epdf* em ) :mpdf ( ) {ep=const_cast<epdf*>( em );};
        void condition ( const vec &cond ) {}
        ~mepdf(){if (owning_ep) delete ep;}
        void from_setting(const Setting &set){
                epdf* e = UI::build<epdf>(set,"epdf");
                ep=     e; 
                owning_ep=true;
        }
};
UIREGISTER(mepdf);

class compositepdf {
protected:
  Array<mpdf*> mpdfs;
  bool owning_mpdfs;
public:
        compositepdf():mpdfs(0){};
        compositepdf(Array<mpdf*> A0, bool own=false){set_elements(A0,own);};
        void set_elements(Array<mpdf*> A0, bool own=false) {mpdfs=A0;owning_mpdfs=own;};
  RV getrv(bool checkoverlap = false);
  void setrvc(const RV &rv, RV &rvc);
  ~compositepdf(){if (owning_mpdfs) for(int i=0;i<mpdfs.length();i++){delete mpdfs(i);}};
};

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


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

  virtual void condition(const vec &val) {it_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 const epdf* _e() 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);
  }
};


}; //namespace
#endif // BDMBASE_H

---