root/library/bdm/base/bdmbase.h @ 610

Revision 610, 33.7 kB (checked in by smidl, 15 years ago)

rv returns name of a scalar at prespecified position

  • Property svn:eol-style set to native
RevLine 
[2]1/*!
[5]2  \file
[508]3  \brief Basic structures of probability calculus: random variables, probability densities, Bayes rule
[5]4  \author Vaclav Smidl.
[2]5
[5]6  -----------------------------------
7  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
8
9  Using IT++ for numerical operations
10  -----------------------------------
11*/
12
[384]13#ifndef BDMBASE_H
14#define BDMBASE_H
[2]15
[351]16#include <map>
[263]17
[190]18#include "../itpp_ext.h"
[357]19#include "../bdmroot.h"
[461]20#include "../shared_ptr.h"
[384]21#include "user_info.h"
[2]22
[340]23using namespace libconfig;
[270]24using namespace itpp;
25using namespace std;
[2]26
[600]27namespace bdm {
[422]28//! Structure of RV, i.e. RVs expanded into a flat list of IDs, used for debugging.
[600]29class str {
30public:
31        //! vector id ids (non-unique!)
32        ivec ids;
33        //! vector of times
34        ivec times;
35        //!Default constructor
36        str ( ivec ids0, ivec times0 ) : ids ( ids0 ), times ( times0 ) {
37                bdm_assert_debug ( times0.length() == ids0.length(), "Incompatible input" );
38        };
[270]39};
[145]40
[270]41/*!
42* \brief Class representing variables, most often random variables
[5]43
[270]44The purpose of this class is to decribe a vector of data. Such description is used for connecting various vectors between each other, see class datalink.
[32]45
[270]46The class is implemented using global variables to assure uniqueness of description:
[2]47
[270]48 In is a vector
49\dot
50digraph datalink {
51rankdir=LR;
52subgraph cluster0 {
53node [shape=record];
[600]54label = "MAP \n std::map<string,int>";
[270]55map [label="{{\"a\"| \"b\" | \"c\"} | {<3> 3 |<1> 1|<2> 2}}"];
56color = "white"
57}
58subgraph cluster1{
59node [shape=record];
[600]60label = "NAMES";
[270]61names [label="{<1> \"b\" | <2> \"c\" | <3>\"a\" }"];
62color = "white"
63}
64subgraph cluster2{
65node [shape=record];
[600]66label = "SIZES";
[270]67labelloc = b;
68sizes [label="{<1>1 |<2> 4 |<3> 1}"];
69color = "white"
70}
71map:1 -> names:1;
72map:1 -> sizes:1;
73map:3 -> names:3;
74map:3 -> sizes:3;
75}
76\enddot
77*/
[32]78
[600]79class RV : public root {
80private:
81        typedef std::map<string, int> str2int_map;
[5]82
[600]83        //! Internal global variable storing sizes of RVs
84        static ivec SIZES;
85        //! Internal global variable storing names of RVs
86        static Array<string> NAMES;
[271]87
[600]88        //! TODO
89        const static int BUFFER_STEP;
[422]90
[600]91        //! TODO
92        static str2int_map MAP;
[422]93
[600]94public:
[422]95
[600]96protected:
97        //! size of the data vector
98        int dsize;
99        //! number of individual rvs
100        int len;
101        //! Vector of unique IDs
102        ivec ids;
103        //! Vector of shifts from current time
104        ivec times;
[532]105
[600]106private:
107        enum enum_dummy {dummy};
[532]108
[600]109        //! auxiliary function used in constructor
110        void init ( const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times );
111        int init ( const string &name, int size );
112        //! Private constructor from IDs, potentially dangerous since all ids must be valid!
113        //! dummy is there to prevent confusion with RV(" string");
114        explicit RV ( const ivec &ids0, enum_dummy dum ) : dsize ( 0 ), len ( ids0.length() ), ids ( ids0 ), times ( zeros_i ( ids0.length() ) ) {
115                dsize = countsize();
116        }
117public:
118        //! \name Constructors
119        //!@{
[271]120
[600]121        //! Full constructor
122        RV ( const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times ) {
123                init ( in_names, in_sizes, in_times );
124        }
[271]125
[600]126        //! Constructor with times=0
127        RV ( const Array<std::string> &in_names, const ivec &in_sizes ) {
128                init ( in_names, in_sizes, zeros_i ( in_names.length() ) );
129        }
[422]130
[600]131        //! Constructor with sizes=1, times=0
132        RV ( const Array<std::string> &in_names ) {
133                init ( in_names, ones_i ( in_names.length() ), zeros_i ( in_names.length() ) );
134        }
[422]135
[600]136        //! Constructor of empty RV
137        RV() : dsize ( 0 ), len ( 0 ), ids ( 0 ), times ( 0 ) {}
[271]138
[600]139        //! Constructor of a single RV with given id
140        RV ( string name, int sz, int tm = 0 );
[271]141
[600]142        // compiler-generated copy constructor is used
143        //!@}
[422]144
[600]145        //! \name Access functions
146        //!@{
[422]147
[600]148        //! State output, e.g. for debugging.
149        friend std::ostream &operator<< ( std::ostream &os, const RV &rv );
[422]150
[604]151        string to_string() {ostringstream o; o << this; return o.str();}
152       
153        //! total size of a random variable
[600]154        int _dsize() const {
155                return dsize;
156        }
[271]157
[600]158        //! access function
159        const ivec& _ids() const {
160                return ids;
161        }
162
163        //! Recount size of the corresponding data vector
164        int countsize() const;
[604]165        //! Vector of cumulative sizes of RV
[600]166        ivec cumsizes() const;
[604]167        //! Number of named parts
[600]168        int length() const {
169                return len;
170        }
171        int id ( int at ) const {
172                return ids ( at );
173        }
174        int size ( int at ) const {
175                return SIZES ( ids ( at ) );
176        }
177        int time ( int at ) const {
178                return times ( at );
179        }
180        std::string name ( int at ) const {
181                return NAMES ( ids ( at ) );
182        }
[610]183        //! returns name of a scalar at position scalat, i.e. it can be in the middle of vector name, in that case it adds "_%d" to it
184        std::string scalarname ( int scalat ) const {
185                bdm_assert_debug(scalat<dsize,"Wrong input index");
186                int id=0;
187                int scalid=0;
188                while (scalid+SIZES(ids(id))<=scalat)  { scalid+=SIZES(ids(id)); id++;};
189                //now id is the id of variable of interest
190                if (size(id)==1)
191                        return  NAMES ( ids ( id ) );
192                else
193                        return  NAMES ( ids ( id ) )+ "_" + num2str(scalat-scalid);
194
195        }
[600]196        void set_time ( int at, int time0 ) {
197                times ( at ) = time0;
198        }
199        //!@}
200
201        //! \name Algebra on Random Variables
202        //!@{
203
204        //! Find indices of self in another rv, \return ivec of the same size as self.
205        ivec findself ( const RV &rv2 ) const;
206        //! Find indices of self in another rv, ignore time, \return ivec of the same size as self.
207        ivec findself_ids ( const RV &rv2 ) const;
208        //! Compare if \c rv2 is identical to this \c RV
209        bool equal ( const RV &rv2 ) const;
210        //! Add (concat) another variable to the current one, \return true if all rv2 were added, false if rv2 is in conflict
211        bool add ( const RV &rv2 );
212        //! Subtract  another variable from the current one
213        RV subt ( const RV &rv2 ) const;
214        //! Select only variables at indices ind
215        RV subselect ( const ivec &ind ) const;
216
217        //! Select only variables at indices ind
218        RV operator() ( const ivec &ind ) const {
219                return subselect ( ind );
220        }
221
222        //! Select from data vector starting at di1 to di2
223        RV operator() ( int di1, int di2 ) const;
224
225        //! Shift \c time by delta.
[604]226        void t_plus ( int delta );
[600]227        //!@}
228
[604]229        //! @{ \name Time manipulation functions
[600]230        //! returns rvs with time set to 0 and removed duplicates
231        RV remove_time() const {
232                return RV ( unique ( ids ), dummy );
233        }
[604]234        //! create new RV from the current one with time shifted by given value
235        RV copy_t(int dt) const {
236                RV tmp=*this;
237                tmp.t_plus(dt);
238                return tmp;
239        }
[600]240        //! return rvs with expanded delayes and sorted in the order of: \f$ [ rv_{0}, rv_{-1}, rv_{
241        RV expand_delayes() const {
242                RV rvt = this->remove_time(); //rv at t=0
243                RV tmp = rvt;
244                int td = mint();
245                for ( int i = -1; i >= td; i-- ) {
[604]246                        rvt.t_plus ( -1 );
[598]247                        tmp.add ( rvt ); //shift u1
248                }
[600]249                return tmp;
250        }
251        //!@}
[271]252
[600]253        //!\name Relation to vectors
254        //!@{
[271]255
[600]256        //! generate \c str from rv, by expanding sizes
257        str tostr() const;
258        //! when this rv is a part of bigger rv, this function returns indices of self in the data vector of the bigger crv.
259        //! Then, data can be copied via: data_of_this = cdata(ind);
260        ivec dataind ( const RV &crv ) const;
261        //! same as dataind but this time crv should not be complete supperset of rv.
262        ivec dataind_part ( const RV &crv ) const;
263        //! generate mutual indices when copying data between self and crv.
264        //! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i)
265        void dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const;
266        //! Minimum time-offset
267        int mint() const {
[603]268                        return times.length()>0 ? min (times) : 0;
[600]269        }
270        //!@}
[357]271
[600]272        /*! \brief UI for class RV (description of data vectors)
[357]273
[600]274        \code
275        rv = {
276            class = "RV"; // class name
277          // UNIQUE IDENTIFIER same names = same variable
278          names = ( "a", "b", "c", ...);   // which will be used e.g. in loggers
[357]279
[600]280          //optional arguments
281          sizes = [1, 2, 3, ...];         // (optional) default = ones()
282          times = [-1, -2, 0, ...];       // time shifts with respect to current time (optional) default = zeros()
283        }
284        \endcode
285        */
286        void from_setting ( const Setting &set );
287
288        // TODO dodelat void to_setting( Setting &set ) const;
289
290        //! Invalidate all named RVs. Use before initializing any RV instances, with care...
291        static void clear_all();
[604]292        //! function for debugging RV related stuff
293        string show_all();
294
[270]295};
[600]296UIREGISTER ( RV );
297SHAREDPTR ( RV );
[32]298
[145]299//! Concat two random variables
[600]300RV concat ( const RV &rv1, const RV &rv2 );
[2]301
[85]302//! Class representing function \f$f(x)\f$ of variable \f$x\f$ represented by \c rv
[600]303class fnc : public root {
304protected:
305        //! Length of the output vector
306        int dimy;
307public:
308        //!default constructor
309        fnc() {};
310        //! function evaluates numerical value of \f$f(x)\f$ at \f$x=\f$ \c cond
311        virtual vec eval ( const vec &cond ) {
312                return vec ( 0 );
313        };
[2]314
[600]315        //! function substitutes given value into an appropriate position
316        virtual void condition ( const vec &val ) {};
[27]317
[600]318        //! access function
319        int dimension() const {
320                return dimy;
321        }
[270]322};
[2]323
[270]324class mpdf;
[7]325
[4]326//! Probability density function with numerical statistics, e.g. posterior density.
[32]327
[600]328class epdf : public root {
329protected:
330        //! dimension of the random variable
331        int dim;
332        //! Description of the random variable
333        RV rv;
[32]334
[600]335public:
336        /*! \name Constructors
337         Construction of each epdf should support two types of constructors:
338        \li empty constructor,
339        \li copy constructor,
[271]340
[600]341        The following constructors should be supported for convenience:
342        \li constructor followed by calling \c set_parameters()
343        \li constructor accepting random variables calling \c set_rv()
[271]344
[600]345         All internal data structures are constructed as empty. Their values (including sizes) will be set by method \c set_parameters(). This way references can be initialized in constructors.
346        @{*/
347        epdf() : dim ( 0 ), rv() {};
348        epdf ( const epdf &e ) : dim ( e.dim ), rv ( e.rv ) {};
349        epdf ( const RV &rv0 ) : dim ( rv0._dsize() ) {
350                set_rv ( rv0 );
351        };
352        void set_parameters ( int dim0 ) {
353                dim = dim0;
354        }
355        //!@}
[271]356
[600]357        //! \name Matematical Operations
358        //!@{
[271]359
[600]360        //! Returns a sample, \f$ x \f$ from density \f$ f_x()\f$
361        virtual vec sample() const {
362                bdm_error ( "not implemented" );
363                return vec();
364        }
[502]365
[600]366        //! Returns N samples, \f$ [x_1 , x_2 , \ldots \ \f$  from density \f$ f_x(rv)\f$
367        virtual mat sample_m ( int N ) const;
[502]368
[600]369        //! Compute log-probability of argument \c val
370        //! In case the argument is out of suport return -Infinity
371        virtual double evallog ( const vec &val ) const {
372                bdm_error ( "not implemented" );
373                return 0.0;
374        }
[502]375
[600]376        //! Compute log-probability of multiple values argument \c val
377        virtual vec evallog_m ( const mat &Val ) const;
[502]378
[600]379        //! Compute log-probability of multiple values argument \c val
380        virtual vec evallog_m ( const Array<vec> &Avec ) const;
[502]381
[600]382        //! Return conditional density on the given RV, the remaining rvs will be in conditioning
383        virtual shared_ptr<mpdf> condition ( const RV &rv ) const;
[271]384
[600]385        //! Return marginal density on the given RV, the remainig rvs are intergrated out
386        virtual shared_ptr<epdf> marginal ( const RV &rv ) const;
[271]387
[600]388        //! return expected value
389        virtual vec mean() const {
390                bdm_error ( "not implemneted" );
391                return vec();
392        }
[271]393
[600]394        //! return expected variance (not covariance!)
395        virtual vec variance() const {
396                bdm_error ( "not implemneted" );
397                return vec();
398        }
[565]399
[600]400        //! Lower and upper bounds of \c percentage % quantile, returns mean-2*sigma as default
401        virtual void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const {
402                vec mea = mean();
403                vec std = sqrt ( variance() );
404                lb = mea - 2 * std;
405                ub = mea + 2 * std;
406        };
407        //!@}
[271]408
[600]409        //! \name Connection to other classes
410        //! Description of the random quantity via attribute \c rv is optional.
411        //! For operations such as sampling \c rv does not need to be set. However, for \c marginalization
412        //! and \c conditioning \c rv has to be set. NB:
413        //! @{
[271]414
[600]415        //!Name its rv
416        void set_rv ( const RV &rv0 ) {
417                rv = rv0;
418        }
[565]419
[600]420        //! True if rv is assigned
421        bool isnamed() const {
422                bool b = ( dim == rv._dsize() );
423                return b;
424        }
[565]425
[600]426        //! Return name (fails when isnamed is false)
427        const RV& _rv() const {
428                bdm_assert_debug ( isnamed(), "" );
429                return rv;
430        }
431        //!@}
[377]432
[600]433        //! \name Access to attributes
434        //! @{
[377]435
[600]436        //! Size of the random variable
437        int dimension() const {
438                return dim;
439        }
440        //! Load from structure with elements:
441        //!  \code
442        //! { rv = {class="RV", names=(...),}; // RV describing meaning of random variable
443        //!   // elements of offsprings
444        //! }
445        //! \endcode
446        //!@}
447        void from_setting ( const Setting &set ) {
448                shared_ptr<RV> r = UI::build<RV> ( set, "rv", UI::optional );
449                if ( r ) {
450                        set_rv ( *r );
[477]451                }
[600]452        }
[377]453
[270]454};
[600]455SHAREDPTR ( epdf );
[32]456
[536]457//! Conditional probability density, e.g. modeling \f$ f( x | y) \f$, where \f$ x \f$ is random variable, \c rv, and \f$ y \f$ is conditioning variable, \c rvc.
[600]458class mpdf : public root {
459protected:
460        //!dimension of the condition
461        int dimc;
462        //! random variable in condition
463        RV rvc;
464private:
465        //! internal epdf, used only as cache to avoid virtual calls of \c _rv() and \c _dimension()
466        epdf* ep;
[461]467
[600]468protected:
469        //! set internal pointer \c ep to point to given \c iepdf
470        void set_ep ( epdf &iepdf ) {
471                ep = &iepdf;
472        }
473        //! set internal pointer \c ep to point to given \c iepdf
474        void set_ep ( epdf *iepdfp ) {
475                ep = iepdfp;
476        }
[271]477
[600]478public:
479        //! \name Constructors
480        //! @{
[2]481
[600]482        mpdf() : dimc ( 0 ), rvc(), ep ( NULL ) { }
[271]483
[600]484        mpdf ( const mpdf &m ) : dimc ( m.dimc ), rvc ( m.rvc ), ep ( m.ep ) { }
485        //!@}
[102]486
[600]487        //! \name Matematical operations
488        //!@{
[461]489
[600]490        //! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv
491        virtual vec samplecond ( const vec &cond ) {
492                bdm_error ( "Not implemented" );
493                return vec();
494        }
[461]495
[600]496        //! Returns \param N samples from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv
497        virtual mat samplecond_m ( const vec &cond, int N );
[461]498
[600]499        //! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
500        virtual double evallogcond ( const vec &dt, const vec &cond ) {
501                bdm_error ( "Not implemented" );
502                return 0.0;
503        }
[201]504
[600]505        //! Matrix version of evallogcond
506        virtual vec evallogcond_m ( const mat &Dt, const vec &cond ) {
507                vec v ( Dt.cols() );
508                for ( int i = 0; i < Dt.cols(); i++ ) {
509                        v ( i ) = evallogcond ( Dt.get_col ( i ), cond );
[489]510                }
[600]511                return v;
512        }
[271]513
[600]514        //! Array<vec> version of evallogcond
515        virtual vec evallogcond_m ( const Array<vec> &Dt, const vec &cond ) {
516                bdm_error ( "Not implemented" );
517                return vec();
518        }
[271]519
[600]520        //! \name Access to attributes
521        //! @{
[461]522
[600]523        const RV& _rv() const {
524                return ep->_rv();
525        }
[604]526        const RV& _rvc() const {
[600]527                return rvc;
528        }
[527]529
[600]530        int dimension() const {
531                return ep->dimension();
532        }
533        int dimensionc() {
534                return dimc;
535        }
[461]536
[600]537        //! Load from structure with elements:
538        //!  \code
539        //! { class = "mpdf_offspring",
540        //!   rv = {class="RV", names=(...),}; // RV describing meaning of random variable
541        //!   rvc= {class="RV", names=(...),}; // RV describing meaning of random variable in condition
542        //!   // elements of offsprings
543        //! }
544        //! \endcode
545        //!@}
546        void from_setting ( const Setting &set );
547        //!@}
[488]548
[600]549        //! \name Connection to other objects
550        //!@{
551        void set_rvc ( const RV &rvc0 ) {
552                rvc = rvc0;
553        }
554        void set_rv ( const RV &rv0 ) {
555                ep->set_rv ( rv0 );
556        }
557        bool isnamed() {
558                return ( ep->isnamed() ) && ( dimc == rvc._dsize() );
559        }
560        //!@}
[270]561};
[600]562SHAREDPTR ( mpdf );
[32]563
[536]564//! Mpdf with internal epdf that is modified by function \c condition
[487]565template <class EPDF>
[600]566class mpdf_internal: public mpdf {
567protected :
568        //! Internal epdf used for sampling
569        EPDF iepdf;
570public:
571        //! constructor
572        mpdf_internal() : mpdf(), iepdf() {
573                set_ep ( iepdf );
574        }
[565]575
[600]576        //! Update \c iepdf so that it represents this mpdf conditioned on \c rvc = cond
577        //! This function provides convenient reimplementation in offsprings
578        virtual void condition ( const vec &cond ) {
579                bdm_error ( "Not implemented" );
580        }
[565]581
[600]582        //!access function to iepdf
583        EPDF& e() {
584                return iepdf;
585        }
[488]586
[600]587        //! Reimplements samplecond using \c condition()
588        vec samplecond ( const vec &cond );
589        //! Reimplements evallogcond using \c condition()
590        double evallogcond ( const vec &val, const vec &cond );
591        //! Efficient version of evallogcond for matrices
592        virtual vec evallogcond_m ( const mat &Dt, const vec &cond );
593        //! Efficient version of evallogcond for Array<vec>
594        virtual vec evallogcond_m ( const Array<vec> &Dt, const vec &cond );
595        //! Efficient version of samplecond
596        virtual mat samplecond_m ( const vec &cond, int N );
[487]597};
598
[270]599/*! \brief DataLink is a connection between two data vectors Up and Down
[2]600
[270]601Up can be longer than Down. Down must be fully present in Up (TODO optional)
602See chart:
603\dot
604digraph datalink {
[377]605  node [shape=record];
606  subgraph cluster0 {
607    label = "Up";
608      up [label="<1>|<2>|<3>|<4>|<5>"];
609    color = "white"
[270]610}
[377]611  subgraph cluster1{
612    label = "Down";
613    labelloc = b;
614      down [label="<1>|<2>|<3>"];
615    color = "white"
[270]616}
617    up:1 -> down:1;
618    up:3 -> down:2;
619    up:5 -> down:3;
620}
621\enddot
[263]622
[270]623*/
[600]624class datalink {
625protected:
626        //! Remember how long val should be
627        int downsize;
[424]628
[600]629        //! Remember how long val of "Up" should be
630        int upsize;
[424]631
[600]632        //! val-to-val link, indices of the upper val
633        ivec v2v_up;
[545]634
[600]635public:
636        //! Constructor
637        datalink() : downsize ( 0 ), upsize ( 0 ) { }
[424]638
[600]639        //! Convenience constructor
640        datalink ( const RV &rv, const RV &rv_up ) {
641                set_connection ( rv, rv_up );
642        }
[271]643
[600]644        //! set connection, rv must be fully present in rv_up
645        void set_connection ( const RV &rv, const RV &rv_up );
[286]646
[600]647        //! set connection using indices
648        void set_connection ( int ds, int us, const ivec &upind );
[424]649
[600]650        //! Get val for myself from val of "Up"
651        vec pushdown ( const vec &val_up ) {
652                vec tmp ( downsize );
653                filldown ( val_up, tmp );
654                return tmp;
655        }
656        //! Get val for vector val_down from val of "Up"
657        void filldown ( const vec &val_up, vec &val_down ) {
658                bdm_assert_debug ( upsize == val_up.length(), "Wrong val_up" );
659                val_down = val_up ( v2v_up );
660        }
661        //! Fill val of "Up" by my pieces
662        void pushup ( vec &val_up, const vec &val ) {
663                bdm_assert_debug ( downsize == val.length(), "Wrong val" );
664                bdm_assert_debug ( upsize == val_up.length(), "Wrong val_up" );
665                set_subvector ( val_up, v2v_up, val );
666        }
667        //! access functions
668        int _upsize() {
669                return upsize;
670        }
671        //! access functions
672        int _downsize() {
673                return downsize;
674        }
[270]675};
[115]676
[598]677/*! Extension of datalink to fill only part of Down
678*/
[600]679class datalink_part : public datalink {
680protected:
[598]681        //! indeces of values in vector downsize
682        ivec v2v_down;
[600]683public:
684        void set_connection ( const RV &rv, const RV &rv_up );
685        //! Get val for vector val_down from val of "Up"
686        void filldown ( const vec &val_up, vec &val_down ) {
687                set_subvector ( val_down, v2v_down, val_up ( v2v_up ) );
688        }
[598]689};
690
691/*! \brief Datalink that buffers delayed values - do not forget to call step()
692
693Up is current data, Down is their subset with possibly delayed values
694*/
[600]695class datalink_buffered: public datalink_part {
696protected:
[598]697        //! History, ordered as \f$[Up_{t-1},Up_{t-2}, \ldots]\f$
698        vec history;
699        //! rv of the history
700        RV Hrv;
[600]701        //! h2v : indeces in down
[598]702        ivec h2v_down;
703        //! h2v : indeces in history
704        ivec h2v_hist;
[603]705        //! v2h: indeces of up too be pushed to h
706        ivec v2h_up;
[600]707public:
708
709        datalink_buffered() : datalink_part(), history ( 0 ), h2v_down ( 0 ), h2v_hist ( 0 ) {};
[598]710        //! push current data to history
[600]711        void step ( const vec &val_up ) {
712                if ( Hrv._dsize() > 0 ) {
713                        history.shift_right ( 0, Hrv._dsize() );
[603]714                  history.set_subvector ( 0, val_up(v2h_up) ); // write U after Drv
[598]715                }
716        }
[600]717        //! Get val for myself from val of "Up"
718        vec pushdown ( const vec &val_up ) {
719                vec tmp ( downsize );
720                filldown ( val_up, tmp );
721                return tmp;
[598]722        }
[600]723
724        void filldown ( const vec &val_up, vec &val_down ) {
725                bdm_assert_debug ( val_down.length() >= downsize, "short val_down" );
726
727                set_subvector ( val_down, v2v_down, val_up ( v2v_up ) ); // copy direct values
728                set_subvector ( val_down, h2v_down, history ( h2v_hist ) ); // copy delayed values
[598]729        }
[600]730
731        void set_connection ( const RV &rv, const RV &rv_up ) {
732                // create link between up and down
733                datalink_part::set_connection ( rv, rv_up );
734
735                // create rvs of history
736                // we can store only what we get in rv_up - everything else is removed
737                ivec valid_ids = rv.findself_ids ( rv_up );
738                RV rv_hist = rv.subselect ( find ( valid_ids >= 0 ) ); // select only rvs that are in rv_up
[603]739          RV rv_hist0 =rv_hist.remove_time();  // these RVs will form history at time =0
740          v2h_up = rv_hist0.dataind(rv_up); // indeces of elements of rv_up to be copied
741          // now we need to know what is needed from Up
[600]742                rv_hist = rv_hist.expand_delayes(); // full regressor - including time 0
[603]743          Hrv=rv_hist.subt(rv_hist0);   // remove time 0
[600]744                history = zeros ( Hrv._dsize() );
745
746                Hrv.dataind ( rv, h2v_hist, h2v_down );
747
748                downsize = v2v_down.length() + h2v_down.length();
749                upsize = v2v_up.length();
750        }
[598]751};
752
753//! buffered datalink from 2 vectors to 1
[600]754class datalink_2to1_buffered {
755protected:
[598]756        datalink_buffered dl1;
757        datalink_buffered dl2;
[600]758public:
[598]759        //! set connection between RVs
[600]760        void set_connection ( const RV &rv, const RV &rv_up1, const RV &rv_up2 ) {
761                dl1.set_connection ( rv, rv_up1 );
762                dl2.set_connection ( rv, rv_up2 );
[598]763        }
[600]764        void filldown ( const vec &val1, const vec &val2, vec &val_down ) {
765                bdm_assert_debug ( val_down.length() >= dl1._downsize() + dl2._downsize(), "short val_down" );
766                dl1.filldown ( val1, val_down );
767                dl2.filldown ( val2, val_down );
[598]768        }
[600]769        void step ( const vec &dt, const vec &ut ) {
770                dl1.step ( dt );
771                dl2.step ( ut );
772        }
[598]773};
774
[600]775
776
[424]777//! Data link with a condition.
[600]778class datalink_m2e: public datalink {
779protected:
780        //! Remember how long cond should be
781        int condsize;
[424]782
[600]783        //!upper_val-to-local_cond link, indices of the upper val
784        ivec v2c_up;
[424]785
[600]786        //!upper_val-to-local_cond link, indices of the local cond
787        ivec v2c_lo;
[192]788
[600]789public:
790        //! Constructor
791        datalink_m2e() : condsize ( 0 ) { }
[545]792
[600]793        //! Set connection between vectors
794        void set_connection ( const RV &rv, const RV &rvc, const RV &rv_up );
[424]795
[600]796        //!Construct condition
797        vec get_cond ( const vec &val_up );
[545]798
[600]799        //! Copy corresponding values to Up.condition
800        void pushup_cond ( vec &val_up, const vec &val, const vec &cond );
[270]801};
[424]802
[192]803//!DataLink is a connection between mpdf and its superordinate (Up)
804//! This class links
[600]805class datalink_m2m: public datalink_m2e {
806protected:
807        //!cond-to-cond link, indices of the upper cond
808        ivec c2c_up;
809        //!cond-to-cond link, indices of the local cond
810        ivec c2c_lo;
[424]811
[600]812public:
813        //! Constructor
814        datalink_m2m() {};
815        //! Set connection between the vectors
816        void set_connection ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) {
817                datalink_m2e::set_connection ( rv, rvc, rv_up );
818                //establish c2c connection
819                rvc.dataind ( rvc_up, c2c_lo, c2c_up );
820                bdm_assert_debug ( c2c_lo.length() + v2c_lo.length() == condsize, "cond is not fully given" );
821        }
[424]822
[600]823        //! Get cond for myself from val and cond of "Up"
824        vec get_cond ( const vec &val_up, const vec &cond_up ) {
825                vec tmp ( condsize );
826                set_subvector ( tmp, v2c_lo, val_up ( v2c_up ) );
827                set_subvector ( tmp, c2c_lo, cond_up ( c2c_up ) );
828                return tmp;
829        }
830        //! Fill
[190]831
[270]832};
[190]833
[270]834/*!
835@brief Class for storing results (and semi-results) of an experiment
[267]836
[270]837This class abstracts logging of results from implementation. This class replaces direct logging of results (e.g. to files or to global variables) by calling methods of a logger. Specializations of this abstract class for specific storage method are designed.
838 */
[600]839class logger : public root {
840protected:
841        //! RVs of all logged variables.
842        Array<RV> entries;
843        //! Names of logged quantities, e.g. names of algorithm variants
844        Array<string> names;
845public:
846        //!Default constructor
847        logger() : entries ( 0 ), names ( 0 ) {}
[267]848
[600]849        //! returns an identifier which will be later needed for calling the \c logit() function
850        //! For empty RV it returns -1, this entry will be ignored by \c logit().
851        virtual int add ( const RV &rv, string prefix = "" ) {
852                int id;
853                if ( rv._dsize() > 0 ) {
854                        id = entries.length();
855                        names = concat ( names, prefix ); // diff
856                        entries.set_length ( id + 1, true );
857                        entries ( id ) = rv;
858                } else {
859                        id = -1;
[477]860                }
[600]861                return id; // identifier of the last entry
862        }
[267]863
[600]864        //! log this vector
865        virtual void logit ( int id, const vec &v ) {
866                bdm_error ( "Not implemented" );
867        };
868        //! log this double
869        virtual void logit ( int id, const double &d ) {
870                bdm_error ( "Not implemented" );
871        };
[267]872
[600]873        //! Shifts storage position for another time step.
874        virtual void step() {
875                bdm_error ( "Not implemneted" );
876        };
[267]877
[600]878        //! Finalize storing information
879        virtual void finalize() {};
[267]880
[600]881        //! Initialize the storage
882        virtual void init() {};
[267]883
[270]884};
[267]885
[270]886/*! \brief Unconditional mpdf, allows using epdf in the role of mpdf.
[190]887
[270]888*/
[600]889class mepdf : public mpdf {
890        //! Internal shared pointer to epdf
891        shared_ptr<epdf> iepdf;
892public:
893        //!Default constructor
894        mepdf() { }
895        //! Set internal shared pointer
896        mepdf ( shared_ptr<epdf> em ) {
897                iepdf = em;
898                set_ep ( *iepdf.get() );
899                dimc = 0;
900        }
[461]901
[600]902        //! empty
903        vec samplecond ( const vec &cond ) {
904                return iepdf->sample();
905        }
906        double evallogcond ( const vec &val, const vec &cond ) {
907                return iepdf->evallog ( val );
908        }
[461]909
[600]910        //! Load from structure with elements:
911        //!  \code
912        //! { class = "mepdf",
913        //!   epdf = {class="epdf_offspring",...}
914        //! }
915        //! \endcode
916        //!@}
917        void from_setting ( const Setting &set );
[270]918};
[600]919UIREGISTER ( mepdf );
920SHAREDPTR ( mepdf );
[115]921
[507]922//! \brief Combines RVs from a list of mpdfs to a single one.
923RV get_composite_rv ( const Array<shared_ptr<mpdf> > &mpdfs, bool checkoverlap = false );
[175]924
[270]925/*! \brief Abstract class for discrete-time sources of data.
[12]926
[603]927The class abstracts operations of:
928\li  data aquisition,
929\li  data-preprocessing, such as  scaling of data,
930\li  data resampling from the task of estimation and control.
[270]931Moreover, for controlled systems, it is able to receive the desired control action and perform it in the next step. (Or as soon as possible).
[12]932
[603]933The DataSource has three main data interaction structures:
934\li input, \f$ u_t \f$,
935\li output \f$ y_t \f$,
936\li data, \f$ d_t=[y_t,u_t, \ldots ]\f$ a collection of all inputs and outputs and possibly some internal variables too.
937
[270]938*/
[32]939
[600]940class DS : public root {
941protected:
[603]942                //! size of data returned by \c getdata()
[600]943        int dtsize;
[603]944                //! size of data
[600]945        int utsize;
[603]946                //!size of output
947                int ytsize;
[600]948        //!Description of data returned by \c getdata().
949        RV Drv;
950        //!Description of data witten by by \c write().
951        RV Urv; //
[603]952                //!Description of output data
953                RV Yrv; //
[600]954        //! Remember its own index in Logger L
955        int L_dt, L_ut;
956public:
957        //! default constructors
[603]958                DS() : Drv(), Urv(),Yrv() {};
[565]959
[609]960        //! Returns maximum number of provided data, by default it is set to maximum allowed length, shorter DS should overload this method! See, MemDS.max_length().
961        virtual int max_length() {return std::numeric_limits< int >::max();}
[600]962        //! Returns full vector of observed data=[output, input]
963        virtual void getdata ( vec &dt ) {
964                bdm_error ( "abstract class" );
965        }
[565]966
[600]967        //! Returns data records at indeces.
968        virtual void getdata ( vec &dt, const ivec &indeces ) {
969                bdm_error ( "abstract class" );
970        }
[565]971
[600]972        //! Accepts action variable and schedule it for application.
973        virtual void write ( vec &ut ) {
974                bdm_error ( "abstract class" );
975        }
[565]976
[600]977        //! Accepts action variables at specific indeces
978        virtual void write ( vec &ut, const ivec &indeces ) {
979                bdm_error ( "abstract class" );
980        }
[32]981
[600]982        //! Moves from \f$ t \f$ to \f$ t+1 \f$, i.e. perfroms the actions and reads response of the system.
983        virtual void step() = 0;
[32]984
[600]985        //! Register DS for logging into logger L
986        virtual void log_add ( logger &L ) {
987                bdm_assert_debug ( dtsize == Drv._dsize(), "invalid DS: dtsize (" + num2str ( dtsize ) + ") different from Drv " + num2str ( Drv._dsize() ) );
988                bdm_assert_debug ( utsize == Urv._dsize(), "invalid DS: utsize (" + num2str ( utsize ) + ") different from Urv " + num2str ( Urv._dsize() ) );
[32]989
[600]990                L_dt = L.add ( Drv, "" );
991                L_ut = L.add ( Urv, "" );
992        }
993        //! Register DS for logging into logger L
994        virtual void logit ( logger &L ) {
995                vec tmp ( Drv._dsize() + Urv._dsize() );
996                getdata ( tmp );
997                // d is first in getdata
998                L.logit ( L_dt, tmp.left ( Drv._dsize() ) );
999                // u follows after d in getdata
1000                L.logit ( L_ut, tmp.mid ( Drv._dsize(), Urv._dsize() ) );
1001        }
1002        //!access function
1003        virtual const RV& _drv() const {
1004                //              return concat (Drv, Urv);// why!!!
1005                return Drv;// why!!!
1006        }
1007        //!access function
1008        const RV& _urv() const {
1009                return Urv;
1010        }
[603]1011                //!access function
1012                const RV& _yrv() const {
1013                        return Yrv;
1014                }
[600]1015        //! set random variables
[603]1016                virtual void set_drv (const  RV &yrv, const RV &urv) {
1017                    Yrv = yrv;
1018                        Drv = concat(yrv,urv);
[600]1019                Urv = urv;
1020        }
[270]1021};
[18]1022
[270]1023/*! \brief Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities.
[32]1024
[283]1025This object represents exact or approximate evaluation of the Bayes rule:
1026\f[
1027f(\theta_t | d_1,\ldots,d_t) = \frac{f(y_t|\theta_t,\cdot) f(\theta_t|d_1,\ldots,d_{t-1})}{f(y_t|d_1,\ldots,d_{t-1})}
1028\f]
1029
1030Access to the resulting posterior density is via function \c posterior().
1031
1032As a "side-effect" it also evaluates log-likelihood of the data, which can be accessed via function _ll().
1033It can also evaluate predictors of future values of \f$y_t\f$, see functions epredictor() and predictor().
1034
1035Alternatively, it can evaluate posterior density conditioned by a known constant, \f$ c_t \f$:
1036\f[
1037f(\theta_t | c_t, d_1,\ldots,d_t) \propto  f(y_t,\theta_t|c_t,\cdot, d_1,\ldots,d_{t-1})
1038\f]
1039
1040The value of \f$ c_t \f$ is set by function condition().
1041
[270]1042*/
[32]1043
[600]1044class BM : public root {
1045protected:
1046        //! Random variable of the data (optional)
1047        RV drv;
1048        //!Logarithm of marginalized data likelihood.
1049        double ll;
1050        //!  If true, the filter will compute likelihood of the data record and store it in \c ll . Set to false if you want to save computational time.
1051        bool evalll;
1052public:
1053        //! \name Constructors
1054        //! @{
[271]1055
[600]1056        BM() : ll ( 0 ), evalll ( true ), LIDs ( 4 ), LFlags ( 4 ) {
1057                LIDs = -1;/*empty IDs*/
1058                LFlags = 0;
1059                LFlags ( 0 ) = 1;  /*log only mean*/
1060        };
1061        BM ( const BM &B ) :  drv ( B.drv ), ll ( B.ll ), evalll ( B.evalll ) {}
1062        //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
1063        //! Prototype: \code BM* _copy_() const {return new BM(*this);} \endcode
1064        virtual BM* _copy_() const {
1065                return NULL;
1066        };
1067        //!@}
[18]1068
[600]1069        //! \name Mathematical operations
1070        //!@{
[271]1071
[600]1072        /*! \brief Incremental Bayes rule
1073        @param dt vector of input data
1074        */
1075        virtual void bayes ( const vec &dt ) = 0;
1076        //! Batch Bayes rule (columns of Dt are observations)
1077        virtual void bayesB ( const mat &Dt );
1078        //! Evaluates predictive log-likelihood of the given data record
1079        //! I.e. marginal likelihood of the data with the posterior integrated out.
1080        virtual double logpred ( const vec &dt ) const {
1081                bdm_error ( "Not implemented" );
1082                return 0.0;
1083        }
[565]1084
[600]1085        //! Matrix version of logpred
1086        vec logpred_m ( const mat &dt ) const {
1087                vec tmp ( dt.cols() );
1088                for ( int i = 0; i < dt.cols(); i++ ) {
1089                        tmp ( i ) = logpred ( dt.get_col ( i ) );
[488]1090                }
[600]1091                return tmp;
1092        }
[32]1093
[600]1094        //!Constructs a predictive density \f$ f(d_{t+1} |d_{t}, \ldots d_{0}) \f$
1095        virtual epdf* epredictor() const {
1096                bdm_error ( "Not implemented" );
1097                return NULL;
1098        };
1099        //!Constructs conditional density of 1-step ahead predictor \f$ f(d_{t+1} |d_{t+h-1}, \ldots d_{t}) \f$
1100        virtual mpdf* predictor() const {
1101                bdm_error ( "Not implemented" );
1102                return NULL;
1103        };
1104        //!@}
[271]1105
[600]1106        //! \name Extension to conditional BM
1107        //! This extension is useful e.g. in Marginalized Particle Filter (\ref bdm::MPF).
1108        //! Alternatively, it can be used for automated connection to DS when the condition is observed
1109        //!@{
[283]1110
[600]1111        //! Name of extension variable
1112        RV rvc;
1113        //! access function
1114        const RV& _rvc() const {
1115                return rvc;
1116        }
[283]1117
[600]1118        //! Substitute \c val for \c rvc.
1119        virtual void condition ( const vec &val ) {
1120                bdm_error ( "Not implemented!" );
1121        }
[283]1122
[600]1123        //!@}
[283]1124
1125
[600]1126        //! \name Access to attributes
1127        //!@{
[271]1128
[600]1129        const RV& _drv() const {
1130                return drv;
1131        }
1132        void set_drv ( const RV &rv ) {
1133                drv = rv;
1134        }
1135        void set_rv ( const RV &rv ) {
1136                const_cast<epdf&> ( posterior() ).set_rv ( rv );
1137        }
1138        double _ll() const {
1139                return ll;
1140        }
1141        void set_evalll ( bool evl0 ) {
1142                evalll = evl0;
1143        }
1144        virtual const epdf& posterior() const = 0;
1145        //!@}
[28]1146
[600]1147        //! \name Logging of results
1148        //!@{
[200]1149
[600]1150        //! Set boolean options from a string, recognized are: "logbounds,logll"
1151        virtual void set_options ( const string &opt ) {
1152                LFlags ( 0 ) = 1;
1153                if ( opt.find ( "logbounds" ) != string::npos ) {
1154                        LFlags ( 1 ) = 1;
1155                        LFlags ( 2 ) = 1;
[477]1156                }
[600]1157                if ( opt.find ( "logll" ) != string::npos ) {
1158                        LFlags ( 3 ) = 1;
1159                }
1160        }
1161        //! IDs of storages in loggers 4:[1=mean,2=lb,3=ub,4=ll]
1162        ivec LIDs;
[190]1163
[600]1164        //! Flags for logging - same size as LIDs, each entry correspond to the same in LIDs
1165        ivec LFlags;
1166        //! Add all logged variables to a logger
1167        virtual void log_add ( logger &L, const string &name = "" ) {
1168                // internal
1169                RV r;
1170                if ( posterior().isnamed() ) {
1171                        r = posterior()._rv();
1172                } else {
1173                        r = RV ( "est", posterior().dimension() );
1174                };
[190]1175
[600]1176                // Add mean value
1177                if ( LFlags ( 0 ) ) LIDs ( 0 ) = L.add ( r, name + "mean_" );
1178                if ( LFlags ( 1 ) ) LIDs ( 1 ) = L.add ( r, name + "lb_" );
1179                if ( LFlags ( 2 ) ) LIDs ( 2 ) = L.add ( r, name + "ub_" );
1180                if ( LFlags ( 3 ) ) LIDs ( 3 ) = L.add ( RV ( "ll", 1 ), name );    //TODO: "local" RV
1181        }
1182        virtual void logit ( logger &L ) {
1183                L.logit ( LIDs ( 0 ), posterior().mean() );
1184                if ( LFlags ( 1 ) || LFlags ( 2 ) ) {  //if one of them is off, its LID==-1 and will not be stored
1185                        vec ub, lb;
1186                        posterior().qbounds ( lb, ub );
1187                        L.logit ( LIDs ( 1 ), lb );
1188                        L.logit ( LIDs ( 2 ), ub );
[477]1189                }
[600]1190                if ( LFlags ( 3 ) ) L.logit ( LIDs ( 3 ), ll );
1191        }
1192        //!@}
1193        void from_setting ( const Setting &set ) {
1194                shared_ptr<RV> r = UI::build<RV> ( set, "drv", UI::optional );
1195                if ( r ) {
1196                        set_drv ( *r );
[488]1197                }
[600]1198                string opt;
1199                if ( !UI::get ( opt, set, "options", UI::optional ) ) {
1200                        set_options ( opt );
[573]1201                }
[600]1202        }
[573]1203
[270]1204};
[32]1205
[529]1206typedef Array<shared_ptr<epdf> > epdf_array;
1207
1208typedef Array<shared_ptr<mpdf> > mpdf_array;
1209
[488]1210template<class EPDF>
[600]1211vec mpdf_internal<EPDF>::samplecond ( const vec &cond ) {
1212        condition ( cond );
[488]1213        vec temp = iepdf.sample();
1214        return temp;
1215}
[339]1216
[488]1217template<class EPDF>
[600]1218mat mpdf_internal<EPDF>::samplecond_m ( const vec &cond, int N ) {
1219        condition ( cond );
1220        mat temp ( dimension(), N );
1221        vec smp ( dimension() );
1222        for ( int i = 0; i < N; i++ ) {
[488]1223                smp = iepdf.sample();
[600]1224                temp.set_col ( i, smp );
[488]1225        }
1226
1227        return temp;
1228}
1229
1230template<class EPDF>
[600]1231double mpdf_internal<EPDF>::evallogcond ( const vec &dt, const vec &cond ) {
[488]1232        double tmp;
[600]1233        condition ( cond );
1234        tmp = iepdf.evallog ( dt );
[488]1235        return tmp;
1236}
1237
1238template<class EPDF>
[600]1239vec mpdf_internal<EPDF>::evallogcond_m ( const mat &Dt, const vec &cond ) {
1240        condition ( cond );
1241        return iepdf.evallog_m ( Dt );
[488]1242}
1243
1244template<class EPDF>
[600]1245vec mpdf_internal<EPDF>::evallogcond_m ( const Array<vec> &Dt, const vec &cond ) {
1246        condition ( cond );
1247        return iepdf.evallog_m ( Dt );
[488]1248}
1249
[254]1250}; //namespace
[384]1251#endif // BDMBASE_H
Note: See TracBrowser for help on using the browser.