- Main Page
- Related Pages
- Classes
- Files

- File List
- File Members

# bdmbase.h

Go to the documentation of this file.

    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                 // compiler-generated copy constructor is used
    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                 // TODO dodelat void to_setting( Setting &set ) const;
    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                 }   //it_assert_debug(isnamed(),""); };
    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                         //establish v2c connection
    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                         //establish c2c connection
    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);   // diff
    00690                                 entries.set_length (id + 1, true);
    00691                                 entries (id) = rv;
    00692                         } else {
    00693                                 id = -1;
    00694                         }
    00695                         return id; // identifier of the last entry
    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                         // d is first in getdata
    00804                         L.logit (L_dt, tmp.left (Drv._dsize()));
    00805                         // u follows after d in getdata
    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;/*empty IDs*/
    00859                         LFlags = 0;
    00860                         LFlags (0) = 1;  /*log only mean*/
    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                         // internal
    00969                         RV r;
    00970                         if (posterior().isnamed()) {
    00971                                 r = posterior()._rv();
    00972                         } else {
    00973                                 r = RV ("est", posterior().dimension());
    00974                         };
    00975
    00976                         // Add mean value
    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);              //TODO: "local" RV
    00981                 }
    00982                 virtual void logit (logger &L) {
    00983                         L.logit (LIDs (0), posterior().mean());
    00984                         if (LFlags (1) || LFlags (2)) {        //if one of them is off, its LID==-1 and will not be stored
    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         // it_assert_debug(std::isfinite(tmp), "Infinite value");
    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 }; //namespace
    01046 #endif // BDMBASE_H