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


  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;

  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;

  static void clear_all();
};
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):dim(rv0._dsize()) {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){
          RV* r = UI::build<RV>(set,"rv");
          if (r){
                  set_rv(*r); 
                  delete r;
          }
  }

};


//TODO Samplecond can be generalized

class mpdf : public root
{
protected:
  int dimc;
  RV rvc;

private:
  shared_ptr<epdf> shep;

public:

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

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


    virtual vec samplecond(const vec &cond);

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

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

    virtual double evallogcond(const vec &dt, const vec &cond);

    virtual vec evallogcond_m(const mat &Dt, const vec &cond);

    virtual vec evallogcond_m(const Array<vec> &Dt, const vec &cond);


    RV _rv() { return shep->_rv(); }
    RV _rvc() { it_assert_debug(isnamed(), ""); return rvc; }
    int dimension() { return shep->dimension(); }
    int dimensionc() { return dimc; }

    epdf *e() { return shep.get(); }

    void set_ep(shared_ptr<epdf> ep) { shep = ep; }

    void from_setting(const Setting &set);

    void set_rvc(const RV &rvc0) { rvc = rvc0; }
    void set_rv(const RV &rv0) { shep->set_rv(rv0); }
    bool isnamed() { return (shep->isnamed()) && (dimc == rvc._dsize()); }
};

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) {
    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():condsize(0) { }

  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 {
public:
    mepdf() { }

    mepdf(shared_ptr<epdf> em) {
        set_ep(em);
        dimc = 0;
    }

    void condition(const vec &cond);

    void from_setting(const Setting &set);
};
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

---