root/bdm/stat/libBM.h @ 176

/*!
  \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
  -----------------------------------
*/

#ifndef BM_H
#define BM_H

#include <itpp/itbase.h>
//#include <std>

using namespace itpp;

//! Structure of RV (used internally)
class str{
public:
        ivec ids;
        ivec times;
        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 {
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 ();

        //! 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;
        //! generate indeces into \param crv data vector that form data vector of self.
        ivec dataind(const RV &crv) const;

        //!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 );


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

class fnc {
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 );
        };

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

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


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

class epdf {
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 sampleN ( int N ) const;
        //! Compute probability of argument \c val
        virtual double eval ( const vec &val ) const {return exp ( this->evalpdflog ( val ) );};

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

        //! Compute log-probability of multiple values argument \c val
        virtual vec evalpdflog_m ( const mat &Val ) const {
                vec x ( Val.cols() );
                for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evalpdflog( Val.get_col(i) ) ;}
                return x;
        }

        //! return expected value
        virtual vec mean() 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 {
protected:
        //! modeled random variable
        RV rv;
        //! random variable in condition
        RV rvc;
        //! pointer to internal epdf
        epdf* ep;
public:

        //! Returns the required moment of the epdf
//      virtual fnc moment ( const int order = 1 );
        //! 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 \param ll is a return value of log-likelihood of the sample.
        virtual vec samplecond (const vec &cond, double &ll ) {this->condition ( cond );
        vec temp= ep->sample();
        ll=ep->evalpdflog ( temp );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 (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->evalpdflog ( smp );}
                return temp;
        };
        //! Update \c ep so that it represents this mpdf conditioned on \c rvc = cond
        virtual void condition ( const vec &cond ) {};

        //! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
        virtual double evalcond ( const vec &dt, const vec &cond ) {this->condition ( cond );return ep->eval ( dt );};

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

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

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

WARNING: the class does not check validity of the \c ep pointer nor its existence.
*/
class mepdf : public mpdf {
public:
        //!Default constructor
        mepdf (epdf &em ) :mpdf ( em._rv(),RV() ) {ep=&em;};
};

//!\brief Abstract composition of pdfs, a base for specific classes 
class compositepdf{
        protected:
                //!Number of mpdfs in the composite
                int n;
                //! Elements of composition
                Array<mpdf*> mpdfs;
                //! Indeces of rvs in common rv
                Array<ivec> rvsinrv;
                //! Indeces of rvc in common rv
                Array<ivec> rvcsinrv;
        public:
                compositepdf(Array<mpdf*> A0): n(A0.length()), mpdfs(A0), rvsinrv(n), rvcsinrv(n){};
                RV getrv(bool checkoverlap=false);
                void setrvc(const RV &rv, RV &rvc);
                void setindices(const RV &rv);
                void setrvcinrv(const RV &rvc, Array<ivec> &rvcind);
};

/*! \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 {
protected:
        //!Observed variables, returned by \c getdata().
        RV Drv;
        //!Action variables, accepted by \c write().
        RV Urv; //
public:
        //! Returns full vector of observed data
        void getdata ( vec &dt );
        //! Returns data records at indeces.
        void getdata ( vec &dt, ivec &indeces );
        //! Accepts action variable and schedule it for application.
        void write ( vec &ut );
        //! Accepts action variables at specific indeces
        void write ( vec &ut, ivec &indeces );
        /*! \brief Method that assigns random variables to the datasource.
        Typically, the datasource will be constructed without knowledge of random variables. This method will associate existing variables with RVs.

        (Inherited from m3k, may be deprecated soon).
        */
        void linkrvs ( RV &drv, RV &urv );

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

};

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

*/

class BM {
protected:
        //!Random variable of the posterior
        RV rv;
        //!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 pointer to the epdf representing posterior density on parameters. Use with care!
        virtual const epdf& _epdf() 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;}
        
        //! Destructor for future use;
        virtual ~BM() {};
        //!access function
        const RV& _rv() const {return 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 {
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;}
};

#endif // BM_H