00001
00013 #ifndef BDMBASE_H
00014 #define BDMBASE_H
00015
00016 #include <map>
00017
00018 #include "../itpp_ext.h"
00019 #include "../bdmroot.h"
00020 #include "../shared_ptr.h"
00021 #include "user_info.h"
00022
00023 using namespace libconfig;
00024 using namespace itpp;
00025 using namespace std;
00026
00027 namespace bdm
00028 {
00029
00030 typedef std::map<string, int> RVmap;
00032 extern ivec RV_SIZES;
00034 extern Array<string> RV_NAMES;
00035
00037 class str
00038 {
00039 public:
00041 ivec ids;
00043 ivec times;
00045 str (ivec ids0, ivec times0) : ids (ids0), times (times0) {
00046 it_assert_debug (times0.length() == ids0.length(), "Incompatible input");
00047 };
00048 };
00049
00088 class RV : public root
00089 {
00090 protected:
00092 int dsize;
00094 int len;
00096 ivec ids;
00098 ivec times;
00099
00100 private:
00102 void init (const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times);
00103 int init (const string &name, int size);
00104 public:
00107
00109 RV (const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times) {
00110 init (in_names, in_sizes, in_times);
00111 }
00112
00114 RV (const Array<std::string> &in_names, const ivec &in_sizes) {
00115 init (in_names, in_sizes, zeros_i (in_names.length()));
00116 }
00117
00119 RV (const Array<std::string> &in_names) {
00120 init (in_names, ones_i (in_names.length()), zeros_i (in_names.length()));
00121 }
00122
00124 RV() : dsize (0), len (0), ids (0), times (0) {}
00125
00127 RV (string name, int sz, int tm = 0);
00128
00129
00131
00134
00136 friend std::ostream &operator<< (std::ostream &os, const RV &rv);
00137
00138 int _dsize() const {
00139 return dsize;
00140 }
00141
00143 int countsize() const;
00144 ivec cumsizes() const;
00145 int length() const {
00146 return len;
00147 }
00148 int id (int at) const {
00149 return ids (at);
00150 }
00151 int size (int at) const {
00152 return RV_SIZES (ids (at));
00153 }
00154 int time (int at) const {
00155 return times (at);
00156 }
00157 std::string name (int at) const {
00158 return RV_NAMES (ids (at));
00159 }
00160 void set_time (int at, int time0) {
00161 times (at) = time0;
00162 }
00164
00167
00169 ivec findself (const RV &rv2) const;
00171 bool equal (const RV &rv2) const;
00173 bool add (const RV &rv2);
00175 RV subt (const RV &rv2) const;
00177 RV subselect (const ivec &ind) const;
00178
00180 RV operator() (const ivec &ind) const {
00181 return subselect (ind);
00182 }
00183
00185 RV operator() (int di1, int di2) const;
00186
00188 void t (int delta);
00190
00193
00195 str tostr() const;
00198 ivec dataind (const RV &crv) const;
00201 void dataind (const RV &rv2, ivec &selfi, ivec &rv2i) const;
00203 int mint() const {
00204 return min (times);
00205 }
00207
00222 void from_setting (const Setting &set);
00223
00224
00226 static void clear_all();
00227 };
00228 UIREGISTER (RV);
00229 SHAREDPTR (RV);
00230
00232 RV concat (const RV &rv1, const RV &rv2);
00233
00235 extern RV RV0;
00236
00238
00239 class fnc : public root
00240 {
00241 protected:
00243 int dimy;
00244 public:
00246 fnc() {};
00248 virtual vec eval (const vec &cond) {
00249 return vec (0);
00250 };
00251
00253 virtual void condition (const vec &val) {};
00254
00256 int dimension() const {
00257 return dimy;
00258 }
00259 };
00260
00261 class mpdf;
00262
00264
00265 class epdf : public root
00266 {
00267 protected:
00269 int dim;
00271 RV rv;
00272
00273 public:
00285 epdf() : dim (0), rv() {};
00286 epdf (const epdf &e) : dim (e.dim), rv (e.rv) {};
00287 epdf (const RV &rv0) : dim (rv0._dsize()) {
00288 set_rv (rv0);
00289 };
00290 void set_parameters (int dim0) {
00291 dim = dim0;
00292 }
00294
00297
00299 virtual vec sample() const {
00300 it_error ("not implemented");
00301 return vec (0);
00302 }
00303
00305 virtual mat sample_m (int N) const;
00306
00309 virtual double evallog (const vec &val) const {
00310 it_error ("not implemented");
00311 return 0.0;
00312 }
00313
00315 virtual vec evallog_m (const mat &Val) const;
00316
00318 virtual vec evallog_m (const Array<vec> &Avec) const;
00319
00321 virtual shared_ptr<mpdf> condition (const RV &rv) const;
00322
00324 virtual shared_ptr<epdf> marginal (const RV &rv) const;
00325
00327 virtual vec mean() const {
00328 it_error ("not implemneted");
00329 return vec (0);
00330 };
00331
00333 virtual vec variance() const {
00334 it_error ("not implemneted");
00335 return vec (0);
00336 };
00338 virtual void qbounds (vec &lb, vec &ub, double percentage = 0.95) const {
00339 vec mea = mean();
00340 vec std = sqrt (variance());
00341 lb = mea - 2 * std;
00342 ub = mea + 2 * std;
00343 };
00345
00351
00353 void set_rv (const RV &rv0) {
00354 rv = rv0;
00355 }
00357 bool isnamed() const {
00358 bool b = (dim == rv._dsize());
00359 return b;
00360 }
00362 const RV& _rv() const {
00363 it_assert_debug (isnamed(), "");
00364 return rv;
00365 }
00367
00370
00372 int dimension() const {
00373 return dim;
00374 }
00382 void from_setting (const Setting &set) {
00383 shared_ptr<RV> r = UI::build<RV> ( set, "rv", UI::optional );
00384 if (r) {
00385 set_rv (*r);
00386 }
00387 }
00388
00389 };
00390 SHAREDPTR(epdf);
00391
00393 class mpdf : public root
00394 {
00395 protected:
00397 int dimc;
00399 RV rvc;
00400 private:
00402 epdf* ep;
00403
00404 protected:
00406 void set_ep (epdf &iepdf) {
00407 ep = &iepdf;
00408 }
00410 void set_ep (epdf *iepdfp) {
00411 ep = iepdfp;
00412 }
00413
00414 public:
00417
00418 mpdf() : dimc (0), rvc(), ep (NULL) { }
00419
00420 mpdf (const mpdf &m) : dimc (m.dimc), rvc (m.rvc), ep (m.ep) { }
00422
00425
00427 virtual vec samplecond (const vec &cond) {it_error ("Not implemented");return vec (0);};
00428
00430 virtual mat samplecond_m (const vec &cond, int N);
00431
00433 virtual double evallogcond (const vec &dt, const vec &cond) {it_error ("Not implemented");return 0.0;}
00434
00436 virtual vec evallogcond_m (const mat &Dt, const vec &cond) {
00437 vec v(Dt.cols());
00438 for(int i=0;i<Dt.cols();i++){v(i)= evallogcond(Dt.get_col(i),cond);}
00439 return v;
00440 }
00441
00443 virtual vec evallogcond_m (const Array<vec> &Dt, const vec &cond) {it_error ("Not implemented");return vec (0);}
00444
00447
00448 RV _rv() const {
00449 return ep->_rv();
00450 }
00451 RV _rvc() {
00452 return rvc;
00453 }
00454
00455 int dimension() const {
00456 return ep->dimension();
00457 }
00458 int dimensionc() {
00459 return dimc;
00460 }
00461
00471 void from_setting (const Setting &set);
00473
00476 void set_rvc (const RV &rvc0) {
00477 rvc = rvc0;
00478 }
00479 void set_rv (const RV &rv0) {
00480 ep->set_rv (rv0);
00481 }
00482 bool isnamed() {
00483 return (ep->isnamed()) && (dimc == rvc._dsize());
00484 }
00486 };
00487 SHAREDPTR(mpdf);
00488
00490 template <class EPDF>
00491 class mpdf_internal: public mpdf
00492 {
00493 protected :
00495 EPDF iepdf;
00496 public:
00498 mpdf_internal() : mpdf(), iepdf() {set_ep (iepdf);}
00501 virtual void condition (const vec &cond) {
00502 it_error ("Not implemented");
00503 };
00505 EPDF& e() {return iepdf;}
00506
00508 vec samplecond (const vec &cond);
00510 double evallogcond (const vec &val, const vec &cond);
00512 virtual vec evallogcond_m (const mat &Dt, const vec &cond);
00514 virtual vec evallogcond_m (const Array<vec> &Dt, const vec &cond);
00516 virtual mat samplecond_m (const vec &cond, int N);
00517 };
00518
00544 class datalink
00545 {
00546 protected:
00548 int downsize;
00549
00551 int upsize;
00552
00554 ivec v2v_up;
00555
00556 public:
00558 datalink() : downsize (0), upsize (0) { }
00560 datalink (const RV &rv, const RV &rv_up) {
00561 set_connection (rv, rv_up);
00562 }
00563
00565 void set_connection (const RV &rv, const RV &rv_up) {
00566 downsize = rv._dsize();
00567 upsize = rv_up._dsize();
00568 v2v_up = rv.dataind (rv_up);
00569
00570 it_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up");
00571 }
00572
00574 void set_connection (int ds, int us, const ivec &upind) {
00575 downsize = ds;
00576 upsize = us;
00577 v2v_up = upind;
00578
00579 it_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up");
00580 }
00581
00583 vec pushdown (const vec &val_up) {
00584 it_assert_debug (upsize == val_up.length(), "Wrong val_up");
00585 return get_vec (val_up, v2v_up);
00586 }
00587
00589 void pushup (vec &val_up, const vec &val) {
00590 it_assert_debug (downsize == val.length(), "Wrong val");
00591 it_assert_debug (upsize == val_up.length(), "Wrong val_up");
00592 set_subvector (val_up, v2v_up, val);
00593 }
00594 };
00595
00597 class datalink_m2e: public datalink
00598 {
00599 protected:
00601 int condsize;
00602
00604 ivec v2c_up;
00605
00607 ivec v2c_lo;
00608
00609 public:
00611 datalink_m2e() : condsize (0) { }
00613 void set_connection (const RV &rv, const RV &rvc, const RV &rv_up) {
00614 datalink::set_connection (rv, rv_up);
00615 condsize = rvc._dsize();
00616
00617 rvc.dataind (rv_up, v2c_lo, v2c_up);
00618 }
00619
00621 vec get_cond (const vec &val_up) {
00622 vec tmp (condsize);
00623 set_subvector (tmp, v2c_lo, val_up (v2c_up));
00624 return tmp;
00625 }
00627 void pushup_cond (vec &val_up, const vec &val, const vec &cond) {
00628 it_assert_debug (downsize == val.length(), "Wrong val");
00629 it_assert_debug (upsize == val_up.length(), "Wrong val_up");
00630 set_subvector (val_up, v2v_up, val);
00631 set_subvector (val_up, v2c_up, cond);
00632 }
00633 };
00634
00637 class datalink_m2m: public datalink_m2e
00638 {
00639 protected:
00641 ivec c2c_up;
00643 ivec c2c_lo;
00644
00645 public:
00647 datalink_m2m() {};
00649 void set_connection (const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up) {
00650 datalink_m2e::set_connection (rv, rvc, rv_up);
00651
00652 rvc.dataind (rvc_up, c2c_lo, c2c_up);
00653 it_assert_debug (c2c_lo.length() + v2c_lo.length() == condsize, "cond is not fully given");
00654 }
00655
00657 vec get_cond (const vec &val_up, const vec &cond_up) {
00658 vec tmp (condsize);
00659 set_subvector (tmp, v2c_lo, val_up (v2c_up));
00660 set_subvector (tmp, c2c_lo, cond_up (c2c_up));
00661 return tmp;
00662 }
00664
00665 };
00666
00672 class logger : public root
00673 {
00674 protected:
00676 Array<RV> entries;
00678 Array<string> names;
00679 public:
00681 logger() : entries (0), names (0) {}
00682
00685 virtual int add (const RV &rv, string prefix = "") {
00686 int id;
00687 if (rv._dsize() > 0) {
00688 id = entries.length();
00689 names = concat (names, prefix);
00690 entries.set_length (id + 1, true);
00691 entries (id) = rv;
00692 } else {
00693 id = -1;
00694 }
00695 return id;
00696 }
00697
00699 virtual void logit (int id, const vec &v) = 0;
00701 virtual void logit (int id, const double &d) = 0;
00702
00704 virtual void step() = 0;
00705
00707 virtual void finalize() {};
00708
00710 virtual void init() {};
00711
00712 };
00713
00717 class mepdf : public mpdf
00718 {
00720 shared_ptr<epdf> iepdf;
00721 public:
00723 mepdf() { }
00725 mepdf (shared_ptr<epdf> em) {
00726 iepdf = em;
00727 set_ep (*iepdf.get());
00728 dimc = 0;
00729 }
00730
00732 vec samplecond(const vec &cond){return iepdf->sample();}
00733 double evallogcond(const vec &val, const vec &cond){return iepdf->evallog(val);}
00734
00742 void from_setting (const Setting &set);
00743 };
00744 UIREGISTER (mepdf);
00745 SHAREDPTR (mepdf);
00746
00748 RV get_composite_rv ( const Array<shared_ptr<mpdf> > &mpdfs, bool checkoverlap = false );
00749
00757 class DS : public root
00758 {
00759 protected:
00760 int dtsize;
00761 int utsize;
00763 RV Drv;
00765 RV Urv;
00767 int L_dt, L_ut;
00768 public:
00770 DS() : Drv(), Urv() {};
00772 virtual void getdata (vec &dt) {
00773 it_error ("abstract class");
00774 };
00776 virtual void getdata (vec &dt, const ivec &indeces) {
00777 it_error ("abstract class");
00778 };
00780 virtual void write (vec &ut) {
00781 it_error ("abstract class");
00782 };
00784 virtual void write (vec &ut, const ivec &indeces) {
00785 it_error ("abstract class");
00786 };
00787
00789 virtual void step() = 0;
00790
00792 virtual void log_add (logger &L) {
00793 it_assert_debug (dtsize == Drv._dsize(), "");
00794 it_assert_debug (utsize == Urv._dsize(), "");
00795
00796 L_dt = L.add (Drv, "");
00797 L_ut = L.add (Urv, "");
00798 }
00800 virtual void logit (logger &L) {
00801 vec tmp (Drv._dsize() + Urv._dsize());
00802 getdata (tmp);
00803
00804 L.logit (L_dt, tmp.left (Drv._dsize()));
00805
00806 L.logit (L_ut, tmp.mid (Drv._dsize(), Urv._dsize()));
00807 }
00809 virtual RV _drv() const {
00810 return concat (Drv, Urv);
00811 }
00813 const RV& _urv() const {
00814 return Urv;
00815 }
00817 virtual void set_drv (const RV &drv, const RV &urv) {
00818 Drv = drv;
00819 Urv = urv;
00820 }
00821 };
00822
00844 class BM : public root
00845 {
00846 protected:
00848 RV drv;
00850 double ll;
00852 bool evalll;
00853 public:
00856
00857 BM() : ll (0), evalll (true), LIDs (4), LFlags (4) {
00858 LIDs = -1;
00859 LFlags = 0;
00860 LFlags (0) = 1;
00861 };
00862 BM (const BM &B) : drv (B.drv), ll (B.ll), evalll (B.evalll) {}
00865 virtual BM* _copy_() const {
00866 return NULL;
00867 };
00869
00872
00876 virtual void bayes (const vec &dt) = 0;
00878 virtual void bayesB (const mat &Dt);
00881 virtual double logpred (const vec &dt) const {
00882 it_error ("Not implemented");
00883 return 0.0;
00884 }
00886 vec logpred_m (const mat &dt) const {
00887 vec tmp (dt.cols());
00888 for (int i = 0; i < dt.cols(); i++) {
00889 tmp (i) = logpred (dt.get_col (i));
00890 }
00891 return tmp;
00892 }
00893
00895 virtual epdf* epredictor() const {
00896 it_error ("Not implemented");
00897 return NULL;
00898 };
00900 virtual mpdf* predictor() const {
00901 it_error ("Not implemented");
00902 return NULL;
00903 };
00905
00910
00912 RV rvc;
00914 const RV& _rvc() const {
00915 return rvc;
00916 }
00917
00919 virtual void condition (const vec &val) {
00920 it_error ("Not implemented!");
00921 };
00922
00924
00925
00928
00929 const RV& _drv() const {
00930 return drv;
00931 }
00932 void set_drv (const RV &rv) {
00933 drv = rv;
00934 }
00935 void set_rv (const RV &rv) {
00936 const_cast<epdf&> (posterior()).set_rv (rv);
00937 }
00938 double _ll() const {
00939 return ll;
00940 }
00941 void set_evalll (bool evl0) {
00942 evalll = evl0;
00943 }
00944 virtual const epdf& posterior() const = 0;
00946
00949
00951 virtual void set_options (const string &opt) {
00952 LFlags (0) = 1;
00953 if (opt.find ("logbounds") != string::npos) {
00954 LFlags (1) = 1;
00955 LFlags (2) = 1;
00956 }
00957 if (opt.find ("logll") != string::npos) {
00958 LFlags (3) = 1;
00959 }
00960 }
00962 ivec LIDs;
00963
00965 ivec LFlags;
00967 virtual void log_add (logger &L, const string &name = "") {
00968
00969 RV r;
00970 if (posterior().isnamed()) {
00971 r = posterior()._rv();
00972 } else {
00973 r = RV ("est", posterior().dimension());
00974 };
00975
00976
00977 if (LFlags (0)) LIDs (0) = L.add (r, name + "mean_");
00978 if (LFlags (1)) LIDs (1) = L.add (r, name + "lb_");
00979 if (LFlags (2)) LIDs (2) = L.add (r, name + "ub_");
00980 if (LFlags (3)) LIDs (3) = L.add (RV ("ll", 1), name);
00981 }
00982 virtual void logit (logger &L) {
00983 L.logit (LIDs (0), posterior().mean());
00984 if (LFlags (1) || LFlags (2)) {
00985 vec ub, lb;
00986 posterior().qbounds (lb, ub);
00987 L.logit (LIDs (1), lb);
00988 L.logit (LIDs (2), ub);
00989 }
00990 if (LFlags (3)) L.logit (LIDs (3), ll);
00991 }
00993 };
00994
00995 typedef Array<shared_ptr<epdf> > epdf_array;
00996
00997 typedef Array<shared_ptr<mpdf> > mpdf_array;
00998
00999 template<class EPDF>
01000 vec mpdf_internal<EPDF>::samplecond (const vec &cond)
01001 {
01002 condition (cond);
01003 vec temp = iepdf.sample();
01004 return temp;
01005 }
01006
01007 template<class EPDF>
01008 mat mpdf_internal<EPDF>::samplecond_m (const vec &cond, int N)
01009 {
01010 condition (cond);
01011 mat temp (dimension(), N);
01012 vec smp (dimension());
01013 for (int i = 0; i < N; i++) {
01014 smp = iepdf.sample();
01015 temp.set_col (i, smp);
01016 }
01017
01018 return temp;
01019 }
01020
01021 template<class EPDF>
01022 double mpdf_internal<EPDF>::evallogcond (const vec &dt, const vec &cond)
01023 {
01024 double tmp;
01025 condition (cond);
01026 tmp = iepdf.evallog (dt);
01027
01028 return tmp;
01029 }
01030
01031 template<class EPDF>
01032 vec mpdf_internal<EPDF>::evallogcond_m (const mat &Dt, const vec &cond)
01033 {
01034 condition (cond);
01035 return iepdf.evallog_m (Dt);
01036 }
01037
01038 template<class EPDF>
01039 vec mpdf_internal<EPDF>::evallogcond_m (const Array<vec> &Dt, const vec &cond)
01040 {
01041 condition (cond);
01042 return iepdf.evallog_m (Dt);
01043 }
01044
01045 };
01046 #endif // BDMBASE_H
