root/bdm/stat/libBM.h @ 267

/*!
  \file
  \brief Bayesian Models (bm) that use Bayes rule to learn from observations
  \author Vaclav Smidl.

  -----------------------------------
  BDM++ - C++ library for Bayesian Decision Making under Uncertainty

  Using IT++ for numerical operations
  -----------------------------------
*/

/*!
\defgroup core Core BDM classes
@{
*/
#ifndef BM_H
#define BM_H


#include "../itpp_ext.h"
//#include <std>

namespace bdm {
        using namespace itpp;
        using std::string;

//! Root class of BDM objects
        class bdmroot {
                virtual void print() {}
        };

//! Structure of RV (used internally), i.e. expanded RVs
        class str {
        public:
                //! vector id ids (non-unique!)
                ivec ids;
                //! vector of times
                ivec times;
                //!Default constructor
                str ( ivec ids0, ivec times0 ) :ids ( ids0 ),times ( times0 ) {
                        it_assert_debug ( times0.length() ==ids0.length(),"Incompatible input" );
                };
        };

        /*!
        * \brief Class representing variables, most often random variables

        * More?...
        */

        class RV :public bdmroot {
        protected:
                //! size = sum of sizes
                int tsize;
                //! len = number of individual rvs
                int len;
                //! Vector of unique IDs
                ivec ids;
                //! Vector of sizes
                ivec sizes;
                //! Vector of shifts from current time
                ivec times;
                //! Array of names
                Array<std::string> names;

        private:
                //! auxiliary function used in constructor
                void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times );
        public:
                //! Full constructor
                RV ( Array<std::string> in_names, ivec in_sizes, ivec in_times );
                //! Constructor with times=0
                RV ( Array<std::string> in_names, ivec in_sizes );
                //! Constructor with sizes=1, times=0
                RV ( Array<std::string> in_names );
                //! Constructor of empty RV
                RV ();
                //! Constructor of a single RV with given id
                RV (string name, int id, int sz=1, int tm=0);

                //! Printing output e.g. for debugging.
                friend std::ostream &operator<< ( std::ostream &os, const RV &rv );

                //! Return number of scalars in the RV.
                int count() const {return tsize;} ;
                //! Return length (number of entries) of the RV.
                int length() const {return len;} ;

                //TODO why not inline and later??

                //! Find indexes of self in another rv, \return ivec of the same size as self.
                ivec findself ( const RV &rv2 ) const;
                //! Compare if \c rv2 is identical to this \c RV
                bool equal ( const RV &rv2 ) const;
                //! Add (concat) another variable to the current one, \return true if all rv2 were added, false if rv2 is in conflict
                bool add ( const RV &rv2 );
                //! Subtract  another variable from the current one
                RV subt ( const RV &rv2 ) const;
                //! Select only variables at indeces ind
                RV subselect ( const ivec &ind ) const;
                //! Select only variables at indeces ind
                RV operator() ( const ivec &ind ) const;
                //! Shift \c time shifted by delta.
                void t ( int delta );
                //! generate \c str from rv, by expanding sizes
                str tostr() const;
                //! when this rv is a part of bigger rv, this function returns indeces of self in the data vector of the bigger crv.
                //! Then, data can be copied via: data_of_this = cdata(ind);
                ivec dataind ( const RV &crv ) const;
                //! generate mutual indeces when copying data betwenn self and crv.
                //! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i)
                void dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const;
                //! Minimum time-offset
                int mint () const {return min ( times );};

                //!access function
                Array<std::string>& _names() {return names;};

                //!access function
                int id ( int at ) {return ids ( at );};
                //!access function
                int size ( int at ) {return sizes ( at );};
                //!access function
                int time ( int at ) {return times ( at );};
                //!access function
                std::string name ( int at ) {return names ( at );};

                //!access function
                void set_id ( int at, int id0 ) {ids ( at ) =id0;};
                //!access function
                void set_size ( int at, int size0 ) {sizes ( at ) =size0; tsize=sum ( sizes );};
                //!access function
                void set_time ( int at, int time0 ) {times ( at ) =time0;};

                //!Assign unused ids to this rv
                void newids();
        };

//! Concat two random variables
        RV concat ( const RV &rv1, const RV &rv2 );

//!Default empty RV that can be used as default argument
        extern RV RV0;

//! Class representing function \f$f(x)\f$ of variable \f$x\f$ represented by \c rv

        class fnc :public bdmroot {
        protected:
                //! Length of the output vector
                int dimy;
        public:
                //!default constructor
                fnc ( int dy ) :dimy ( dy ) {};
                //! function evaluates numerical value of \f$f(x)\f$ at \f$x=\f$ \c cond
                virtual vec eval ( const vec &cond ) {
                        return vec ( 0 );
                };

                //! function substitutes given value into an appropriate position
                virtual void condition ( const vec &val ) {};

                //! access function
                int _dimy() const{return dimy;}

                //! Destructor for future use;
                virtual ~fnc() {};
        };

        class mpdf;

//! Probability density function with numerical statistics, e.g. posterior density.

        class epdf :public bdmroot {
        protected:
                //! Identified of the random variable
                RV rv;
        public:
                //!default constructor
                epdf() :rv ( ) {};

                //!default constructor
                epdf ( const RV &rv0 ) :rv ( rv0 ) {};

//      //! Returns the required moment of the epdf
//      virtual vec moment ( const int order = 1 );

                //! Returns a sample, \f$x\f$ from density \f$epdf(rv)\f$
                virtual vec sample () const =0;
                //! Returns N samples from density \f$epdf(rv)\f$
                virtual mat sample_m ( int N ) const;

                //! Compute log-probability of argument \c val
                virtual double evallog ( const vec &val ) const =0;

                //! Compute log-probability of multiple values argument \c val
                virtual vec evallog_m ( const mat &Val ) const {
                        vec x ( Val.cols() );
                        for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog ( Val.get_col ( i ) ) ;}
                        return x;
                }
                //! Return conditional density on the given RV, the remaining rvs will be in conditioning
                virtual mpdf* condition ( const RV &rv ) const  {it_warning ( "Not implemented" ); return NULL;}
                //! Return marginal density on the given RV, the remainig rvs are intergrated out
                virtual epdf* marginal ( const RV &rv ) const {it_warning ( "Not implemented" ); return NULL;}

                //! return expected value
                virtual vec mean() const =0;

                //! return expected variance (not covariance!)
                virtual vec variance() const = 0;

                //! Destructor for future use;
                virtual ~epdf() {};
                //! access function, possibly dangerous!
                const RV& _rv() const {return rv;}
                //! modifier function - useful when copying epdfs
                void _renewrv ( const RV &in_rv ) {rv=in_rv;}
                //!
        };


//! Conditional probability density, e.g. modeling some dependencies.
//TODO Samplecond can be generalized

        class mpdf : public bdmroot {
        protected:
                //! modeled random variable
                RV rv;
                //! random variable in condition
                RV rvc;
                //! pointer to internal epdf
                epdf* ep;
        public:

                //! 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 
                virtual vec samplecond ( const vec &cond) {
                        this->condition ( cond );
                        vec temp= ep->sample();
                        return temp;
                };
                //! 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 \param ll is a return value of log-likelihood of the sample.
                virtual mat samplecond_m ( const vec &cond, vec &ll, int N ) {
                        this->condition ( cond );
                        mat temp ( rv.count(),N ); vec smp ( rv.count() );
                        for ( int i=0;i<N;i++ ) {smp=ep->sample() ;temp.set_col ( i, smp );ll ( i ) =ep->evallog ( smp );}
                        return temp;
                };
                //! Update \c ep so that it represents this mpdf conditioned on \c rvc = cond
                virtual void condition ( const vec &cond ) {it_error ( "Not implemented" );};

                //! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
                virtual double evallogcond ( const vec &dt, const vec &cond ) {
                        double tmp; this->condition ( cond );tmp = ep->evallog ( dt );          it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); return tmp;
                };

                //! Matrix version of evallogcond
                virtual vec evallogcond_m ( const mat &Dt, const vec &cond ) {this->condition ( cond );return ep->evallog_m ( Dt );};

                //! Destructor for future use;
                virtual ~mpdf() {};

                //! Default constructor
                mpdf ( const RV &rv0, const RV &rvc0 ) :rv ( rv0 ),rvc ( rvc0 ) {};
                //! access function
                RV _rvc() const {return rvc;}
                //! access function
                RV _rv() const {return rv;}
                //!access function
                epdf& _epdf() {return *ep;}
                //!access function
                epdf* _e() {return ep;}
        };

        /*! \brief DataLink is a connection between two data vectors Up and Down

        Up can be longer than Down. Down must be fully present in Up (TODO optional)
        See chart:
        \dot
        digraph datalink {
                node [shape=record];
                subgraph cluster0 {
                        label = "Up";
                up [label="<1>|<2>|<3>|<4>|<5>"];
                        color = "white"
        }
                subgraph cluster1{
                        label = "Down";
                        labelloc = b;
                down [label="<1>|<2>|<3>"];
                        color = "white"
        }
            up:1 -> down:1;
            up:3 -> down:2;
            up:5 -> down:3;
        }
        \enddot

        */
        class datalink_e2e {
        protected:
                //! Remember how long val should be
                int valsize;
                //! Remember how long val of "Up" should be
                int valupsize;
                //! val-to-val link, indeces of the upper val
                ivec v2v_up;
        public:
                //! Constructor
                datalink_e2e ( const RV &rv, const RV &rv_up ) :
                                valsize ( rv.count() ), valupsize ( rv_up.count() ), v2v_up ( rv.dataind ( rv_up ) )  {
                        it_assert_debug ( v2v_up.length() ==valsize,"rv is not fully in rv_up" );
                }
                //! Get val for myself from val of "Up"
                vec get_val ( const vec &val_up ) {it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" ); return get_vec ( val_up,v2v_up );}
                //! Fill val of "Up" by my pieces
                void fill_val ( vec &val_up, const vec &val ) {
                        it_assert_debug ( valsize==val.length(),"Wrong val" );
                        it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" );
                        set_subvector ( val_up, v2v_up, val );
                }
        };

//! data link between
        class datalink_m2e: public datalink_e2e {
        protected:
                //! Remember how long cond should be
                int condsize;
                //!upper_val-to-local_cond link, indeces of the upper val
                ivec v2c_up;
                //!upper_val-to-local_cond link, ideces of the local cond
                ivec v2c_lo;

        public:
                //! Constructor
                datalink_m2e ( const RV &rv,  const RV &rvc, const RV &rv_up ) :
                                datalink_e2e ( rv,rv_up ), condsize ( rvc.count() ) {
                        //establish v2c connection
                        rvc.dataind ( rv_up, v2c_lo, v2c_up );
                }
                //!Construct condition
                vec get_cond ( const vec &val_up ) {
                        vec tmp ( condsize );
                        set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) );
                        return tmp;
                }
                void fill_val_cond ( vec &val_up, const vec &val, const vec &cond ) {
                        it_assert_debug ( valsize==val.length(),"Wrong val" );
                        it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" );
                        set_subvector ( val_up, v2v_up, val );
                        set_subvector ( val_up, v2c_up, cond );
                }
        };
//!DataLink is a connection between mpdf and its superordinate (Up)
//! This class links
        class datalink_m2m: public datalink_m2e {
        protected:
                //!cond-to-cond link, indeces of the upper cond
                ivec c2c_up;
                //!cond-to-cond link, indeces of the local cond
                ivec c2c_lo;
        public:
                //! Constructor
                datalink_m2m ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) :
                                datalink_m2e ( rv, rvc, rv_up ) {
                        //establish c2c connection
                        rvc.dataind ( rvc_up, c2c_lo, c2c_up );
                        it_assert_debug ( c2c_lo.length() +v2c_lo.length() ==condsize, "cond is not fully given" );
                }
                //! Get cond for myself from val and cond of "Up"
                vec get_cond ( const vec &val_up, const vec &cond_up ) {
                        vec tmp ( condsize );
                        set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) );
                        set_subvector ( tmp,c2c_lo,cond_up ( c2c_up ) );
                        return tmp;
                }
                //! Fill

        };

        /*!
        @brief Class for storing results (and semi-results) of an experiment

        This 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.
         */
        class logger : public bdmroot {
                protected:
                //! RVs of all logged variables.
                        Array<RV> entries;
                //! Names of logged quantities, e.g. names of algorithm variants
                        Array<string> names;
                public:
                //!Default constructor
                        logger ( ) : entries ( 0 ),names ( 0 ) {}

                //! returns an identifier which will be later needed for calling the log() function
                        virtual int add ( const RV &rv, string name="" ) {
                                int id=entries.length();
                                names=concat ( names, name ); // diff
                                entries.set_length ( id+1,true );
                                entries ( id ) = rv;
                                return id; // identifier of the last entry
                        }

                //! log this vector
                        virtual void logit ( int id, const vec &v ) =0;

                //! Shifts storage position for another time step.
                        virtual void step() =0;

                //! Finalize storing information
                        virtual void finalize() {};

                //! Initialize the storage
                        virtual void init() {};

                //! for future use
                        virtual ~logger() {};
        };

        /*! \brief Unconditional mpdf, allows using epdf in the role of mpdf.

        */
        class mepdf : public mpdf {
        public:
                //!Default constructor
                mepdf ( const epdf* em ) :mpdf ( em->_rv(),RV() ) {ep=const_cast<epdf*> ( em );};
                void condition ( const vec &cond ) {}
        };

//!\brief Abstract composition of pdfs, will be used for specific classes
//!this abstract class is common to epdf and mpdf
        class compositepdf {
        protected:
                //!Number of mpdfs in the composite
                int n;
                //! Elements of composition
                Array<mpdf*> mpdfs;
        public:
                compositepdf ( Array<mpdf*> A0 ) : n ( A0.length() ), mpdfs ( A0 ) {};
                //! find common rv, flag \param checkoverlap modifies whether overlaps are acceptable
                RV getrv ( bool checkoverlap=false );
                //! common rvc of all mpdfs is written to rvc
                void setrvc ( const RV &rv, RV &rvc );
        };

        /*! \brief Abstract class for discrete-time sources of data.

        The class abstracts operations of: (i) data aquisition, (ii) data-preprocessing, (iii) scaling of data, and (iv) data resampling from the task of estimation and control.
        Moreover, for controlled systems, it is able to receive the desired control action and perform it in the next step. (Or as soon as possible).

        */

        class DS : public bdmroot {
        protected:
                //!Observed variables, returned by \c getdata().
                RV Drv;
                //!Action variables, accepted by \c write().
                RV Urv; //
                //! Remember its own index in Logger L
                int L_dt, L_ut;
        public:
                DS() :Drv ( RV0 ),Urv ( RV0 ) {};
                DS ( const RV &Drv0, const RV &Urv0 ) :Drv ( Drv0 ),Urv ( Urv0 ) {};
                //! Returns full vector of observed data=[output, input]
                virtual void getdata ( vec &dt ) {it_error ( "abstract class" );};
                //! Returns data records at indeces.
                virtual void getdata ( vec &dt, const ivec &indeces ) {it_error ( "abstract class" );};
                //! Accepts action variable and schedule it for application.
                virtual void write ( vec &ut ) {it_error ( "abstract class" );};
                //! Accepts action variables at specific indeces
                virtual void write ( vec &ut, const ivec &indeces ) {it_error ( "abstract class" );};

                //! Moves from \f$t\f$ to \f$t+1\f$, i.e. perfroms the actions and reads response of the system.
                virtual void step() =0;

                //! Register DS for logging into logger L
                virtual void log_add ( logger &L ) {
                        L_dt=L.add ( Drv,"" );
                        L_ut=L.add ( Urv,"" );
                }
                //! Register DS for logging into logger L
                virtual void logit ( logger &L ) {
                        vec tmp(Drv.count()+Urv.count());
                        getdata(tmp);
                        // d is first in getdata
                        L.logit ( L_dt,tmp.left ( Drv.count() ) );
                        // u follows after d in getdata
                        L.logit ( L_ut,tmp.mid ( Drv.count(), Urv.count() ) );
                }
                //!access function 
                virtual RV _drv() const {return concat(Drv,Urv);}
                //!access function 
                const RV& _urv() const {return Urv;}
        };

        /*! \brief Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities.

        */

        class BM :public bdmroot {
        protected:
                //!Random variable of the posterior
                RV rv;
                //! Random variable of the data (optional)
                RV drv;
                //!Logarithm of marginalized data likelihood.
                double ll;
                //!  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.
                bool evalll;
        public:

                //!Default constructor
                BM ( const RV &rv0, double ll0=0,bool evalll0=true ) :rv ( rv0 ), ll ( ll0 ),evalll ( evalll0 ) {//Fixme: test rv
                };
                //!Copy constructor
                BM ( const BM &B ) : rv ( B.rv ), ll ( B.ll ), evalll ( B.evalll ) {}

                /*! \brief Incremental Bayes rule
                @param dt vector of input data
                */
                virtual void bayes ( const vec &dt ) = 0;
                //! Batch Bayes rule (columns of Dt are observations)
                virtual void bayesB ( const mat &Dt );
                //! Returns a reference to the epdf representing posterior density on parameters.
                virtual const epdf& _epdf() const =0;

                //! Returns a pointer to the epdf representing posterior density on parameters. Use with care!
                virtual const epdf* _e() const =0;

                //! Evaluates predictive log-likelihood of the given data record
                //! I.e. marginal likelihood of the data with the posterior integrated out.
                virtual double logpred ( const vec &dt ) const{it_error ( "Not implemented" );return 0.0;}
                //! Matrix version of logpred
                vec logpred_m ( const mat &dt ) const{vec tmp ( dt.cols() );for ( int i=0;i<dt.cols();i++ ) {tmp ( i ) =logpred ( dt.get_col ( i ) );}return tmp;}

                //!Constructs a predictive density (marginal density on data)
                virtual epdf* predictor ( const RV &rv ) const {it_error ( "Not implemented" );return NULL;};

                //! Destructor for future use;
                virtual ~BM() {};
                //!access function
                const RV& _rv() const {return rv;}
                //!access function 
                const RV& _drv() const {return drv;}
                //!set drv
                void set_drv(const RV &rv){drv=rv;}
                //!access function
                double _ll() const {return ll;}
                //!access function
                void set_evalll ( bool evl0 ) {evalll=evl0;}

                //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
                //! Prototype: BM* _copy_(){BM Tmp*=new Tmp(this*);  return Tmp; }
                virtual BM* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
        };

        /*!
        \brief Conditional Bayesian Filter

        Evaluates conditional filtering density \f$f(rv|rvc,data)\f$ for a given \c rvc which is specified in each step by calling function \c condition.

        This is an interface class used to assure that certain BM has operation \c condition .

        */

        class BMcond :public bdmroot {
        protected:
                //! Identificator of the conditioning variable
                RV rvc;
        public:
                //! Substitute \c val for \c rvc.
                virtual void condition ( const vec &val ) =0;
                //! Default constructor
                BMcond ( RV &rv0 ) :rvc ( rv0 ) {};
                //! Destructor for future use
                virtual ~BMcond() {};
                //! access function
                const RV& _rvc() const {return rvc;}
        };

}; //namespace
/*! @} */
#endif // BM_H