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

Revision 604, 33.0 kB (checked in by smidl, 15 years ago)

change of syntax of RV

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